In [1]:
import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import sys
import cPickle
import matplotlib.units
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
from common_functions import *
import matplotlib_helper_functions as mplhf
import driven_optical_matter_functions as domf
%matplotlib inline
#%config InlineBackend.figure_format = 'retina'
#import seaborn as sns
In [2]:
os.chdir("C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Data")
In [3]:
#matplotlib.rcParams
In [4]:
matplotlib.rcParams['font.family']
matplotlib.rcParams['mathtext.fontset']
Out[4]:
In [5]:
matplotlib.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
In [6]:
#matplotlib.rcParams['font.family'] = 'Vera'
matplotlib.rcParams['figure.dpi'] = 200.0
#matplotlib.rcParams['axes.labelsize'] = 15
#matplotlib.rcParams['axes.titlesize'] = 18
#matplotlib.rcParams['xtick.labelsize'] = 'large'
#matplotlib.rcParams['ytick.labelsize'] = 'large'
matplotlib.rcParams['mathtext.fontset'] = 'stixsans'
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = 'stixsans'
matplotlib.rcParams['legend.handletextpad'] = 0.1
matplotlib.rcParams['legend.numpoints'] = 1
matplotlib.rcParams['legend.scatterpoints'] = 1
matplotlib.rcParams['legend.borderaxespad'] = 0.9
#matplotlib.rcParams['mathtext.default'] = 'regular'
#matplotlib.rcParams['font.family'] = 'serif'
#matplotlib.rcParams['font.serif'] = 'Times'
#matplotlib.rcParams['text.usetex'] = True
#matplotlib.rcParams['text.latex.unicode'] = True
matplotlib.rc('text.latex', preamble=[r'\usepackage{berasans}',
r'\renewcommand*\familydefault{\sfdefault}',
r'\usepackage[T1]{fontenc}',
r'\usepackage{amsmath}'])
In [7]:
# For Matplotlib 2.0
# matplotlib.rcParams['xtick.direction'] = 'in'
# matplotlib.rcParams['ytick.direction'] = 'in'
# matplotlib.rcParams['xtick.top'] = True
# matplotlib.rcParams['ytick.right'] = True
In [8]:
def add_axes_label_inches(ax, (right, down), string, **kwargs):
"""A helper function to add a text lable a specific distance from the top left corner in inches
This function was created so that when constructing multi-panel figures for a manuscript all the axes labels
(e.g. a, b, c, d...) are a consistent distance from the corner of the axis. The relative coordinates used for text
labels can result in a label being a different physical distance compared to subplots of different sizes.
:param fig: A mpl.figure object that contains the axis you want to add text to
:param ax: A mpl.axes object that you want to add text to
:param (right,down): Distance in inches you want the text label to be from the top left corner of the axis
:param string: The text string you want displayed at that location."""
fig = ax.get_figure()
fig_size = fig.get_size_inches()
ax_bbox = ax.get_position()
ax_rect_inches = ax_bbox.x0*fig_size[0], ax_bbox.y0*fig_size[1], ax_bbox.x1*fig_size[0], ax_bbox.y1*fig_size[1]
text_location_inches = [right, ax_rect_inches[3]-ax_rect_inches[1]-down]
text_position_rel_coors = text_location_inches[0]/(ax_rect_inches[2]-ax_rect_inches[0]), text_location_inches[1]/(ax_rect_inches[3]-ax_rect_inches[1])
return ax.text(text_position_rel_coors[0], text_position_rel_coors[1], string, transform=ax.transAxes, va='top', ha='left', **kwargs)
In [9]:
fig = plt.figure(figsize=[3,4])
ax = plt.subplot()
add_axes_label_inches(ax, (0.3,0.3), 'a', fontsize=14)
plt.close()
In [10]:
import IPython.display
IPython.display.Image(filename="C:/Users/Scherer Lab E/Box Sync/Ring/Figure Data and Scripts/Figures/fig 1.PNG")
Out[10]:
Schematics of experimental setup and optical trap over (a) the coverslip glass and (b) a Au nanoplate. (c) Image of two Ag nanoparticles trapped over glass. (d) Image of two Ag nanoparticles trapped over a Au nanoplate. The curved arrows in (c) and (d) indicate the direction of rotation of the nanoparticles.
In [11]:
f = open('Fig_2.pkl', 'r')
notebook_info, all_082214_data_gl, all_082214_data_pl, um_conv = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [12]:
%matplotlib inline
from matplotlib.colors import LogNorm
# Make custom colormaps to fit the colors for over glass and over plate used in
# the rest of the paper. Always goes from white to the colors selected below
red_color=np.array([186,62,66])
blue_color=np.array([68,102,158])
custom_colorbar_red=matplotlib.colors.LinearSegmentedColormap.from_list('custom_colorbar_red', [(1,1,1),red_color/256.0])
custom_colorbar_blue=matplotlib.colors.LinearSegmentedColormap.from_list('custom_colorbar_blue', [(1,1,1),blue_color/256.0])
um_conv=6.5/60/1.6/2
y_label = [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$\frac{3\pi}{2}$", r"$2\pi$"]
particle_pos_and_polar_transform_fig=plt.figure(figsize=[3,3])
ax1=plt.subplot(2,2,1)
im1=plt.hist2d(all_082214_data_gl['d_x'].values*um_conv,all_082214_data_gl['d_y'].values*um_conv, bins=100, range=[[-5.5,5.5],[-5.5,5.5]], norm=LogNorm(),
cmap=custom_colorbar_blue)
im1 = im1[3]
plt.axis('equal')
plt.xlabel(u'X (µm)')
plt.ylabel(u'Y (µm)')
x_y_ticks=[-5,0,5]
plt.xticks(x_y_ticks,x_y_ticks)
plt.yticks(x_y_ticks,x_y_ticks)
plt.xlim([-6,6])
plt.ylim([-6,6])
ax1.xaxis.labelpad = 0
ax1.yaxis.labelpad = 0
ax2=plt.subplot(2,2,2, sharey=ax1)
im2 = plt.hist2d(all_082214_data_pl['d_x'].values*um_conv,
all_082214_data_pl['d_y'].values*um_conv, bins=100,
range=[[-5.5,5.5],[-5.5,5.5]], norm=LogNorm(),
cmap=custom_colorbar_red)
im2 = im2[3]
clim_2 = im2.get_clim()
print clim_2
plt.setp(ax2.get_yticklabels(), visible=False)
plt.axis('equal')
plt.xlabel(u'X (µm)')
plt.xticks(x_y_ticks,x_y_ticks)
plt.yticks(x_y_ticks,x_y_ticks)
#plt.yticks(ax2.get_xticklabels(),ax2.get_xticks()[:-1])
plt.xlim([-6,6])
plt.ylim([-6,6])
ax2.xaxis.labelpad = 0
ax3=plt.subplot(2,2,3)
# Shift theta so that it follows the conventional unit circle
ax3_theta = (all_082214_data_gl['theta'].values*np.pi/180)+np.pi/2
ax3_theta = ax3_theta%(2*np.pi)
im3 = plt.hist2d(all_082214_data_gl['r'].values*um_conv,
ax3_theta, bins=100, range=[[3,6],[0,2*np.pi]],
norm=LogNorm(), cmap=custom_colorbar_blue)
im3 = im3[3]
im3.set_clim(1, clim_2[1])
y_label = [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$\frac{3\pi}{2}$", r"$2\pi$"]
plt.xlabel(u'Radius (µm)')
plt.ylabel(r'Radians')
plt.xlim([2.5,6.5])
polar_x_ticks=[3,4,5,6]
plt.xticks(polar_x_ticks, polar_x_ticks)
plt.yticks( np.array([0,0.5,1,1.5,2])*np.pi,y_label)
ax3.xaxis.labelpad = 0
ax3.yaxis.labelpad = 0
#plt.yticks([0,90,180,270,360])
ax3.tick_params(axis='y', which='major', labelsize=12)
#plt.colorbar()
ax4=plt.subplot(2,2,4, sharey=ax3)
# Shift theta so that it follows the conventional unit circle
ax4_theta = (all_082214_data_pl['theta'].values*np.pi/180)+np.pi/2
ax4_theta = ax4_theta%(2*np.pi)
im4=plt.hist2d(all_082214_data_pl['r'].values*um_conv,
ax4_theta, range=[[4,6],[0,2*np.pi]],
bins=100, norm=LogNorm(),cmap=custom_colorbar_red)
im4 = im4[3]
clim_4 = im4.get_clim()
# im4.set_clim(1, clim_1[1])
plt.xlabel(u'Radius (µm)')
plt.xlim([2.5,6.5])
plt.ylim([0,2*np.pi])
plt.xticks(polar_x_ticks, polar_x_ticks)
plt.yticks( np.array([0,0.5,1,1.5,2])*np.pi,y_label)
plt.setp(ax4.get_yticklabels(), visible=False)
ax4.xaxis.labelpad = 0
#plt.colorbar()
im1.set_clim(1, clim_4[1])
im2.set_clim(1, clim_4[1])
im3.set_clim(1, clim_4[1])
#plt.pcolor()
#particle_pos_and_polar_transform_fig.colorbar()
#print ax1.bbox
#print particle_pos_and_polar_transform_fig.get_figheight()
cbar_ax = particle_pos_and_polar_transform_fig.add_axes([1, 0.14, 0.05, .8])
cbar=particle_pos_and_polar_transform_fig.colorbar(im1,cax=cbar_ax)
cbar_ax = particle_pos_and_polar_transform_fig.add_axes([1.17, 0.14, 0.05, .8])
cbar=particle_pos_and_polar_transform_fig.colorbar(im4,cax=cbar_ax, use_gridspec=True)
cbar.set_label(r'Probability Density')
particle_pos_and_polar_transform_fig.tight_layout()
particle_pos_and_polar_transform_fig.subplots_adjust(wspace=0.00,hspace=0.45)
particle_pos_and_polar_transform_fig=plt.gcf()
# ax1.text(0.07, .80, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
# ax3.text(0.07, .80, 'b', fontsize=14, weight='extra bold', transform=ax3.transAxes)
add_axes_label_inches(ax1, (0.075,0.075), r'a', fontsize=12, weight='extra bold')
add_axes_label_inches(ax3, (0.075,0.075), r'b', fontsize=12, weight='extra bold')
plt.close()
plt.show()
Simulation Trajectory of L=2¶
This data is a sample trajectory of L=2 from Nishant's simulation. From Nishant:
field_intensity is an array of the field intensity. x and y are the x and y axes for plotting the field intensity. x1, y1 and x2, y2 are the positions of a nanoparticle for intensity I_0 and 4*I_0, respectively. The trajectories are for L=2.
In [13]:
import scipy.io as sio
import re
traj_sample_dat = sio.loadmat('trajectories_L2.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_traj_sample_dat = {}
for key, val in traj_sample_dat.iteritems():
if key in skip_list:
continue
elif key in ['x','y']:
new_traj_sample_dat[key] = val[0, :]
elif key in ['field_intensity']:
new_traj_sample_dat[key] = val
else:
new_traj_sample_dat[key] = val[:, 0]
#print key,val.shape
traj_sample_dat = new_traj_sample_dat
In [14]:
def add_four_polar_coords(ax):
"""A helper function that adds 0pi, pi/2, pi, and 3pi/2 axis labels
on the inside of an axis following the unit circle.
:param ax: A mpl.axes object you want to add polar labels to
"""
# Add the polar labels
ax.text(0.95, 0.5, '$\mathbf{\mathregular{0}}$', fontdict=dict(color='white', size='large'),
weight='extra bold', transform=ax.transAxes, va='center', ha='right')
ax.text(0.5, 0.95, '$\mathbf{\pi/\mathregular{2}}$', fontdict=dict(color='white', size='large'),
weight='extra bold', transform=ax.transAxes, va='top', ha='center')
ax.text(0.05, 0.5, '$\mathbf{\pi}$', fontdict=dict(color='white', size='large'),
weight='extra bold', transform=ax.transAxes, va='center', ha='left')
ax.text(0.5, 0.05, r'$\mathbf{\mathregular{3}\pi/\mathregular{2}}$', fontdict=dict(color='white', size='large'),
weight='extra bold', transform=ax.transAxes, va='bottom', ha='center')
# Change ticks in graph
xticklines = ax.get_xticklines()
yticklines = ax.get_yticklines()
xcenter = len(xticklines)/2
ycenter = len(yticklines)/2
xticklines[xcenter].set(color='white', markeredgewidth=1.5, markersize=4, zorder=0)
xticklines[xcenter-1].set(color='white', markeredgewidth=1.5, markersize=4)
yticklines[ycenter].set(color='white', markeredgewidth=1.5, markersize=4)
yticklines[ycenter-1].set(color='white', markeredgewidth=1.5, markersize=4)
In [15]:
force_fig = plt.figure(figsize=(2,4))
# Plot electric field intensity
field_x = traj_sample_dat['x']/10**3
field_y = traj_sample_dat['y']/10**3
field_int = traj_sample_dat['field_intensity']
ax5 = plt.subplot(211)
ax6 = plt.subplot(212)
ax5.pcolormesh(field_x, field_y, field_int, cmap='viridis')
ax5.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
ax5.set_ylabel(u'Position (µm)')
ax5.set_xlabel(u'Position (µm)')
ax6.pcolormesh(field_x, field_y, field_int, cmap='viridis')
ax6.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
ax6.set_ylabel(u'Position (µm)')
ax6.set_xlabel(u'Position (µm)')
# Plot the sample trajectory
x_pos = traj_sample_dat['x1']/10**3
y_pos = traj_sample_dat['y1']/10**3
ax5.plot(x_pos, y_pos, c="#C44E52")
text_c = ax5.text(0.03, .87, 'c', fontsize=14, weight='extra bold', transform=ax5.transAxes)
text_c.set_bbox(dict(boxstyle="square,pad=0.1", color='white', alpha=0.7))
x_pos = traj_sample_dat['x2']/10**3
y_pos = traj_sample_dat['y2']/10**3
ax6.plot(x_pos, y_pos, c="#C44E52")
text_d = ax6.text(0.03, .87, 'd', fontsize=14, weight='extra bold', transform=ax6.transAxes)
text_d.set_bbox(dict(boxstyle="square,pad=0.1", color='white', alpha=0.7))
ax5.set_aspect('equal')
ax6.set_aspect('equal')
plt.subplots_adjust(hspace=0.45)
plt.close()
Combine the Figures into multipanel¶
In [16]:
from matplotlib.colors import LogNorm
from matplotlib.gridspec import GridSpec
# Make custom colormaps to fit the colors for over glass and over plate used in
# the rest of the paper. Always goes from white to the colors selected below
red_color=np.array([186,62,66])
blue_color=np.array([68,102,158])
custom_colorbar_red=matplotlib.colors.LinearSegmentedColormap.from_list('custom_colorbar_red', [(1,1,1),red_color/256.0])
custom_colorbar_blue=matplotlib.colors.LinearSegmentedColormap.from_list('custom_colorbar_blue', [(1,1,1),blue_color/256.0])
um_conv=6.5/60/1.6/2
y_label = [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$\frac{3\pi}{2}$", r"$2\pi$"]
particle_pos_and_polar_transform_fig=plt.figure(figsize=[3.375,5.54455])
gs1 = GridSpec(2,2, wspace=0, hspace=0.35)
gs1.update(top=1.0, bottom=0.58, right=0.86)
gs_cbar = GridSpec(1,2, wspace=0.3)
gs_cbar.update(left=0.89, right=0.95, bottom=0.58, top=1.0)
ax1=plt.subplot(gs1[0,0])
im1=plt.hist2d(all_082214_data_gl['d_x'].values*um_conv,all_082214_data_gl['d_y'].values*um_conv, bins=100, range=[[-5.5,5.5],[-5.5,5.5]], norm=LogNorm(),
cmap=custom_colorbar_blue, rasterized=True)
im1 = im1[3]
plt.axis('equal')
plt.xlabel(u'x (µm)')
plt.ylabel(u'y (µm)')
x_y_ticks=[-5,0,5]
plt.xticks(x_y_ticks,x_y_ticks)
plt.yticks(x_y_ticks,x_y_ticks)
plt.minorticks_on()
plt.xlim([-6,6])
plt.ylim([-6,6])
ax1.xaxis.labelpad = 0
ax1.yaxis.labelpad = 0
ax2=plt.subplot(gs1[0,1], sharey=ax1)
im2 = plt.hist2d(all_082214_data_pl['d_x'].values*um_conv,
all_082214_data_pl['d_y'].values*um_conv, bins=100,
range=[[-5.5,5.5],[-5.5,5.5]], norm=LogNorm(),
cmap=custom_colorbar_red, rasterized=True)
im2 = im2[3]
plt.setp(ax2.get_yticklabels(), visible=False)
plt.axis('equal')
plt.xlabel(u'x (µm)')
plt.xticks(x_y_ticks,x_y_ticks)
plt.yticks(x_y_ticks,x_y_ticks)
plt.minorticks_on()
#plt.yticks(ax2.get_xticklabels(),ax2.get_xticks()[:-1])
plt.xlim([-6,6])
plt.ylim([-6,6])
ax2.xaxis.labelpad = 0
ax3=plt.subplot(gs1[1,0])
ax3_theta = (all_082214_data_gl['theta'].values*np.pi/180)+np.pi/2
ax3_theta = ax3_theta%(2*np.pi)
# The data was extending beyond the the axes borders on the top. To
# fix this I am going to trim off the data for the last bin near 2pi
ax3_theta_trimmed_index = ax3_theta < 1.95*np.pi
im3 = plt.hist2d(all_082214_data_gl['r'].values[ax3_theta_trimmed_index]*um_conv,
ax3_theta[ax3_theta_trimmed_index], bins=100, range=[[3,6],[0,2*np.pi]],
norm=LogNorm(), cmap=custom_colorbar_blue, rasterized=True)
im3 = im3[3]
y_label = [r"0", r"$\frac{\pi}{\mathregular{2}}$", r"$\pi$", r"$\frac{\mathregular{3}\pi}{\mathregular{2}}$", r"2$\pi$"]
plt.xlabel(u'Radius (µm)')
plt.ylabel(r'Radians')
plt.xlim([2.5,6.5])
polar_x_ticks=[3,4,5,6]
plt.xticks(polar_x_ticks, polar_x_ticks)
plt.yticks(np.array([0,0.5,1,1.5,2])*np.pi,y_label)
ax3.xaxis.labelpad = 0
ax3.yaxis.labelpad = 0
#plt.yticks([0,90,180,270,360])
ax3.tick_params(axis='y', which='major', labelsize=12)
#plt.colorbar()
ax4=plt.subplot(gs1[1,1], sharey=ax3)
ax4_theta = (all_082214_data_pl['theta'].values*np.pi/180)+np.pi/2
ax4_theta = ax4_theta%(2*np.pi)
# The data was extending beyond the the axes borders on the top. To
# fix this I am going to trim off the data for the last bin near 2pi
ax4_theta_trimmed_index = ax4_theta < 1.95*np.pi
im4=plt.hist2d(all_082214_data_pl['r'].values[ax4_theta_trimmed_index]*um_conv,
ax4_theta[ax4_theta_trimmed_index], range=[[4,6],[0,2*np.pi]],
bins=100, norm=LogNorm(),cmap=custom_colorbar_red, rasterized=True)
im4 = im4[3]
clim_4 = im4.get_clim()
plt.xlabel(u'Radius (µm)')
plt.xlim([2.5,6.5])
plt.ylim([0,2*np.pi])
plt.xticks(polar_x_ticks, polar_x_ticks)
plt.yticks( np.array([0,0.5,1,1.5,2])*np.pi,y_label)
plt.setp(ax4.get_yticklabels(), visible=False)
ax4.xaxis.labelpad = 0
#plt.colorbar()
im1.set_clim(1, clim_4[1])
im2.set_clim(1, clim_4[1])
im3.set_clim(1, clim_4[1])
# ax1.text(0.03, .87, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
# ax3.text(0.03, .87, 'b', fontsize=14, weight='extra bold', transform=ax3.transAxes)
add_axes_label_inches(ax1, (0.09,0.068), '(a)', fontsize=10)
add_axes_label_inches(ax3, (0.09,0.068), '(b)', fontsize=10)
#plt.pcolor()
#particle_pos_and_polar_transform_fig.colorbar()
#print ax1.bbox
#print particle_pos_and_polar_transform_fig.get_figheight()
cbar_ax_gl = plt.subplot(gs_cbar[0])
cbar=particle_pos_and_polar_transform_fig.colorbar(im1,cax=cbar_ax_gl)
plt.setp(cbar.ax.yaxis.get_ticklabels(), visible=False)
cbar_ax_pl = plt.subplot(gs_cbar[1])
cbar=particle_pos_and_polar_transform_fig.colorbar(im4,cax=cbar_ax_pl, use_gridspec=False)
cbar.set_ticks([10**0, 10**1, 10**2])
cbar.set_ticklabels([r'$\mathregular{10}^\mathregular{0}$', r'$\mathregular{10}^\mathregular{1}$', r'$\mathregular{10}^\mathregular{2}$'])
cbar.set_label(r'Probability Density', labelpad=1.7)
#particle_pos_and_polar_transform_fig.tight_layout()
#particle_pos_and_polar_transform_fig.subplots_adjust(wspace=0.00,hspace=0.45)
#particle_pos_and_polar_transform_fig=plt.gcf()
gs2 = GridSpec(1,2)
gs2.update(top=0.58, bottom=0.10, wspace=0.06, right=1.0)
ax5 = plt.subplot(gs2[0])
# Plot electric field intensity
field_x = traj_sample_dat['x']/10**3
field_y = traj_sample_dat['y']/10**3
field_int = traj_sample_dat['field_intensity']
ax5 = plt.subplot(gs2[0])
ax6 = plt.subplot(gs2[1])
ax5.pcolormesh(field_x, field_y, field_int, cmap='viridis', rasterized=True)
ax5.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
ax5.set_ylabel(u'y (µm)')
ax5.set_xlabel(u'x (µm)')
ax5.set_yticks([-1,-0.5,0,0.5,1])
ax5.set_yticklabels(['-1','-0.5','0','0.5','1'])
ax5.set_xticks([-1,-0.5,0,0.5,1])
ax5.set_xticklabels(['-1','-0.5','0','0.5','1'])
add_four_polar_coords(ax5)
ax6.pcolormesh(field_x, field_y, field_int, cmap='viridis', rasterized=True)
ax6.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
#ax6.set_ylabel(r'y (µm)')
ax6.set_xlabel(u'x (µm)')
ax6.set_xticks([-1,-0.5,0,0.5,1])
ax6.set_xticklabels(['-1','-0.5','0','0.5','1'])
ax6.set_yticks([-1,-0.5,0,0.5,1])
plt.setp(ax6.get_ymajorticklabels(), visible=False)
#ax6.set_yticklabels(['-1','-0.5','0','0.5','1'])
add_four_polar_coords(ax6)
# Plot the sample trajectory
x_pos = traj_sample_dat['x1']/10**3
y_pos = traj_sample_dat['y1']/10**3
ax5.plot(x_pos, y_pos, c="k", rasterized=True)
text_c = add_axes_label_inches(ax5, (0.09,0.068), '(c)', fontsize=10, color='white')
# text_c.set_bbox(dict(boxstyle="square,pad=0.1", color='white'))
x_pos = traj_sample_dat['x2']/10**3
y_pos = traj_sample_dat['y2']/10**3
ax6.plot(x_pos, y_pos, c="k", rasterized=True)
text_d = add_axes_label_inches(ax6, (0.09,0.068), '(d)', fontsize=10, color='white')
# text_d.set_bbox(dict(boxstyle="square,pad=0.1", color='white'))
ax5.set_aspect('equal')
ax6.set_aspect('equal')
particle_pos_and_polar_transform_fig.set_size_inches(3.34646,4.65)
#top=0.58, bottom=0.10, wspace=0.06, right=1.0
gs1.tight_layout(particle_pos_and_polar_transform_fig, rect=[0, 0.4038, 0.78, 0.998], w_pad=0, h_pad=0.35, pad=0)
gs_cbar.tight_layout(particle_pos_and_polar_transform_fig, rect=[0.80, 0.5, 1, 1], w_pad=0.3, h_pad=0.35, pad=0)
gs2.tight_layout(particle_pos_and_polar_transform_fig, rect=[0, 0, 0.99, 0.462], pad=0, w_pad=0.6)
# Get the Colorbar axes to line up with their graphs
bbox_bot_ax = ax4.get_position()
bbox_top_ax = ax1.get_position()
for ax in [cbar_ax_gl, cbar_ax_pl]:
cur_bbox = ax.get_position()
cur_bbox.y0 = bbox_bot_ax.y0
cur_bbox.y1 = bbox_top_ax.y1
ax.set_position(cur_bbox)
plt.show()
In [17]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
particle_pos_and_polar_transform_fig.savefig(save_dir+"\Fig_2.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
particle_pos_and_polar_transform_fig.savefig(save_dir+"\Fig_2.pdf", dpi=800)
Probability densities of particles in the ring traps over glass (blue) or over plate(red). (a) Distribution of the particle positions in cartesian coordinates in the lab frame. (b) Distributions of particle positions in polar coordinates. All the distributions are log normalized.
In [18]:
f = open('Fig_3.pkl', 'r')
notebook_info, frame_rate, sp_gl_plot, sp_pl_plot, sp_gl, sp_pl, TAMSD_gl_plot, TAMSD_pl_plot, lin_reg_params_glass, xf_gl, yf_gl, lin_reg_params_plate, xf_pl, yf_pl, TAMSD_gl_1_x, TAMSD_gl, TAMSD_gl_1_prams, TAMSD_pl_1_x, TAMSD_pl, TAMSD_pl_1_prams, frame_rate = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [19]:
'''Define functions to help calculate the anglular trajectory'''
def angular_traj(data_frame,x_cent,y_cent):
grp_tracks=data_frame.groupby('track id')
ang_trajs=[]
for name,grp in grp_tracks:
single_traj_angle=[0]
for pos in range(len(grp)-1):
angle_diff=calc_angle(grp['x pos'].iloc[pos+1],grp['y pos'].iloc[pos+1],x_cent,y_cent)-calc_angle(grp['x pos'].iloc[pos],grp['y pos'].iloc[pos],x_cent,y_cent)
if angle_diff>180:
angle_diff=angle_diff-360
if angle_diff<-180:
angle_diff=angle_diff+360
single_traj_angle.append(angle_diff)
ang_trajs.append(single_traj_angle)
return ang_trajs
In [20]:
sp_gl_plot=sp_gl[5:]
sp_pl_plot=sp_pl[5:]
ang_trajs_and_tamsd=plt.figure(figsize=[4.5,2.5])
# ang_trajs_and_tamsd=plt.figure(figsize=[10,8])
ax1=plt.subplot(121)
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
y_ticks = np.array([0,5,10,15,20])
y_label = [r"$0$", r"$5\pi$", r"$10\pi$", r"$15\pi$", r"$20\pi$"]
for num,L in enumerate(sp_gl_plot):
time = np.arange(len(L))/frame_rate
angle = np.cumsum(angular_traj(L,xf_gl,yf_gl))*(np.pi/180)
fit_angle = lin_reg_params_glass[num+5][0]*(np.arange(len(L))/frame_rate)
plt.plot(time, angle, 'o', label=labels[num], color=markercolors[num], markersize=1.5, markeredgewidth=0)
plt.plot(time, fit_angle,'-', color='k', lw=.9, zorder=3)
plt.legend(loc='upper left', labelspacing=-0.5, frameon=False, prop={'size':7})
plt.xlabel('Time (s)')
plt.ylabel('Radians')
plt.yticks(y_ticks*np.pi,y_label)
plt.title('Over Glass')
ax1.xaxis.labelpad=0
ax1.yaxis.labelpad=0
#plt.setp(ax1.get_xticklabels(), visible=False)
ax2=plt.subplot(122)
y_ticks = np.array([0,30,60,90,120])
y_label = [r"$0$", r"$30\pi$", r"$60\pi$", r"$90\pi$", r"$120\pi$"]
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
for num,L in enumerate(sp_pl_plot):
time = np.arange(len(L))/frame_rate
angle = np.cumsum(angular_traj(L,xf_pl,yf_pl))*(np.pi/180)
fit_angle = lin_reg_params_plate[num+5][0]*(np.arange(len(L))/frame_rate)
plt.plot(time, angle, 'o', label=labels[num], color=markercolors[num], markersize=1.5, markeredgewidth=0)
plt.plot(time, fit_angle,'-', color='k', lw=.9, zorder=3)
#plt.ylim([-0.5*np.pi,28*np.pi])
#plt.xlim([0,1.5])
#plt.legend(loc='upper left', labelspacing=0.1, frameon=False, prop={'size':7})
#plt.title('Single Particle Trajectories over Plate')
plt.xlabel('Time (s)')
plt.ylabel('Radians')
plt.title('Over Au Plate')
plt.yticks(y_ticks*np.pi,y_label)
ax2.xaxis.labelpad=0
ax2.yaxis.labelpad=0
ax1.text(0.03, .87, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
ax2.text(0.03, .87, 'b', fontsize=14, weight='extra bold', transform=ax2.transAxes)
plt.tight_layout()
#plt.subplots_adjust(wspace= 0.4, hspace= 0.3)
#sns.despine()
plt.close()
plt.show()
Dynamic information about single Ag nanoparticles in the optical vortex trap for various substrates and angular drive l. Top panels show the trajectory of single Ag nanoparticle rotation with time over glass (a) and over the Au nanoplate (b). Bottom panels show the time averaged mean squared displacement of the single particle rotating over glass (c) and over the Au nanoplate (d). The dashed black curves in (a) and (b) show linear fits of the rotation versus time while the dashed black curves in (c) and (d) show fits to the TAMSD using eq (1).
Angular Speed and TAMSD Simulations¶
This is from Nishant's 2 particle simulations of the ring trap.
In [21]:
import scipy.io as sio
import re
ang_sim_dat = sio.loadmat('angular_displacement_msd.mat')
ang_msd_fit_dat = sio.loadmat('msd_fitting_params.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_sim_dat
new_ang_sim_dat = {}
for key, val in ang_sim_dat.iteritems():
if key in skip_list:
continue
new_ang_sim_dat[key] = val[:,0]
ang_sim_dat = new_ang_sim_dat
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_ang_msd_fit_dat = {}
for key, val in ang_msd_fit_dat.iteritems():
if key in skip_list:
continue
new_ang_msd_fit_dat[key] = val[0, :]
ang_msd_fit_dat = new_ang_msd_fit_dat
Fit the Simulation Data¶
In [22]:
import scipy
# Fit trajectory over glass
sim_lin_reg_params_glass=[]
for L in range(6):
x = ang_sim_dat['time1']*10**4
y = ang_sim_dat['tang_l'+str(L)+'e1']
func=lambda x, *params: params*x
params,cov = scipy.optimize.curve_fit(func,x,y,p0=[1])
sim_lin_reg_params_glass.append(params)
#print sim_lin_reg_params_glass
# Fit trajectory over plate
sim_lin_reg_params_plate=[]
for L in range(6):
x = ang_sim_dat['time1']*10**4
y = ang_sim_dat['tang_l'+str(L)+'e2']
func=lambda x, *params: params*x
params,cov = scipy.optimize.curve_fit(func,x,y,p0=[1])
sim_lin_reg_params_plate.append(params)
In [23]:
import fnmatch
def plot_sim_results(axis, ang_sim_dat, time_key, regex):
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
keys = sorted(ang_sim_dat.keys())
match_keys = fnmatch.filter(keys, regex)
for num, key in enumerate(match_keys):
axis.plot(ang_sim_dat[time_key]*10**4, ang_sim_dat[key], 'o',
c=markercolors[num], markersize=1.7, markeredgewidth=0, rasterized=True)
axis.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
axis.xaxis.labelpad=0
axis.yaxis.labelpad=0
def plot_sim_fits(axis, ang_sim_dat, fit_params, time_key, msd=False):
x = ang_sim_dat[time_key]*10**4
if msd == False:
fit_function = lambda x,*prams: x*prams[0]
real_fit_params = fit_params
for L in range(6):
axis.plot(x, fit_function(x, real_fit_params[L]), '-k', rasterized=True, lw=.9, zorder=3)
elif msd == True:
fit_function=lambda x,*prams: prams[0]*x**2+prams[1]*x
real_fit_params = []
for v in fit_params:
real_fit_params.append(v)
for L in range(6):
axis.plot(x, fit_function(x, real_fit_params[L][0], real_fit_params[L][1]), '-k', rasterized=True, lw=.9, zorder=3)
In [24]:
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
plt.figure(figsize=[5,4])
ax3 = plt.subplot(221)
plot_sim_results(ax3, ang_sim_dat, 'time1', 'tang_l*e1')
plot_sim_fits(ax3, ang_sim_dat, sim_lin_reg_params_glass, 'time1', msd=False)
# plt.plot(time, angle, 'o', label=labels[num], color=markercolors[num], markersize=2, markeredgewidth=0)
# plt.plot(time, fit_angle,'-', color='k', lw=.9, zorder=3)
plt.legend(labels, loc='upper left', labelspacing=0.0, frameon=False, prop={'size':7})
plt.xlabel('Time ($10^{-4}$s)')
plt.ylabel('Radians')
plt.title('Intensity = I$_{o}$')
plt.xlim(0,6.5)
y_ticks = np.array([0,1,2,3])
y_label = [r"$0$", r"$1\pi$", r"$2\pi$", r"$3\pi$"]
plt.yticks(y_ticks*np.pi,y_label)
ax4 = plt.subplot(222)
plot_sim_results(ax4, ang_sim_dat, 'time1', 'tang_l*e2')
plot_sim_fits(ax4, ang_sim_dat, sim_lin_reg_params_plate, 'time1', msd=False)
plt.xlabel('Time ($10^{-4}$s)')
plt.ylabel('Radians')
plt.title('Intensity = 4I$_{o}$')
plt.xlim(0,6.5)
y_ticks = np.array([0,5,10])
y_label = [r"$0$", r"$5\pi$", r"$10\pi$"]
plt.yticks(y_ticks*np.pi,y_label)
ax3.text(0.03, .87, 'a', fontsize=14, weight='extra bold', transform=ax3.transAxes)
ax4.text(0.03, .87, 'b', fontsize=14, weight='extra bold', transform=ax4.transAxes)
plt.tight_layout()
plt.close()
Combine All Plots¶
In [25]:
matplotlib.rcParams['xtick.labelsize'] = 'small'
matplotlib.rcParams['ytick.labelsize'] = 'small'
title_fontsize = matplotlib.rcParams['axes.titlesize']
legend_fontsize = matplotlib.rcParams['legend.fontsize']
# Plot Experiment Trajectories
sp_gl_plot=sp_gl[5:]
sp_pl_plot=sp_pl[5:]
ang_trajs=plt.figure(figsize=[3.34646, 2.677168])
ax1=plt.subplot(221)
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
y_ticks = np.array([0,5,10,15,20])*np.pi
y_label = [r"0", r"5$\pi$", r"10$\pi$", r"15$\pi$", r"20$\pi$"]
for num,L in enumerate(sp_gl_plot):
time = np.arange(len(L))/frame_rate
angle = np.cumsum(angular_traj(L,xf_gl,yf_gl))*(np.pi/180)
fit_angle = lin_reg_params_glass[num+5][0]*(np.arange(len(L))/frame_rate)
plt.plot(time, angle, 'o', rasterized=True, label=labels[num], color=markercolors[num], markersize=1.7, markeredgewidth=0)
plt.plot(time, fit_angle,'-', rasterized=True, color='k', lw=.9, zorder=3)
# plt.legend(loc='upper left', labelspacing=0,
# frameon=False, prop={'size':7}, ncol=2, columnspacing=0.2)
# #plt.xlabel('Time (s)')
plt.ylabel('Radians')
plt.yticks(y_ticks,y_label)
#plt.title('Over Glass')
ax1.text(0.53, .95, r'Over Glass', transform=ax1.transAxes,
va='top', ha='center', fontsize=7.5, usetex=True)
ax1.xaxis.labelpad=0
ax1.yaxis.labelpad=0
ylim_1 = ax1.get_ylim()
ax1.set_ylim(ymax=1.1*ylim_1[1])
plt.setp(ax1.get_xmajorticklabels(), visible=False)
ax1.text(10, -0.2*np.pi, labels[0], transform=ax1.transData, va='top', ha='left', fontsize=6.5, usetex=True)
ax1.text(10, 3.15*np.pi, labels[1], transform=ax1.transData, va='top', ha='left', fontsize=6.5, usetex=True)
ax1.text(10, 7.4*np.pi, labels[2], transform=ax1.transData, va='top', ha='left', fontsize=6.5, usetex=True)
ax1.text(10, 9.70*np.pi, labels[3], transform=ax1.transData, va='top', ha='left', fontsize=6.5, usetex=True)
ax1.text(9.2, 14.7*np.pi, labels[4], transform=ax1.transData, va='top', ha='left', fontsize=6.5, usetex=True)
ax1.text(9.5, 17.9*np.pi, labels[5], transform=ax1.transData, va='top', ha='left', fontsize=6.5, usetex=True)
#plt.setp(ax1.get_xticklabels(), visible=False)
ax3=plt.subplot(223)
y_ticks = np.array([0,30,60,90,120])*np.pi
y_label = [r"0", r"30$\pi$", r"60$\pi$", r"90$\pi$", r"120$\pi$"]
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
for num,L in enumerate(sp_pl_plot):
time = np.arange(len(L))/frame_rate
angle = np.cumsum(angular_traj(L,xf_pl,yf_pl))*(np.pi/180)
fit_angle = lin_reg_params_plate[num+5][0]*(np.arange(len(L))/frame_rate)
plt.plot(time, angle, 'o', rasterized=True, label=labels[num], color=markercolors[num], markersize=1.7, markeredgewidth=0)
plt.plot(time, fit_angle,'-', rasterized=True, color='k', lw=.9, zorder=3)
#plt.ylim([-0.5*np.pi,28*np.pi])
#plt.xlim([0,1.5])
#plt.legend(loc='upper left', labelspacing=0.1, frameon=False, prop={'size':7})
#plt.title('Single Particle Trajectories over Plate')
plt.xlabel('Time (s)')
plt.ylabel('Radians')
#plt.title('Over Au Plate')
ax3.text(0.53, .95, r'Over Au Plate', transform=ax3.transAxes,
va='top', ha='center', fontsize=7.5, usetex=True)
plt.yticks(y_ticks,y_label)
#ax3.set_ylim(ymax=24*np.pi)
ax3.xaxis.labelpad=0
ax3.yaxis.labelpad=0
ylim_3 = ax3.get_ylim()
ax3.set_ylim(ymax=1.0*ylim_3[1])
# ax1.text(0.03, .87, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
# ax2.text(0.03, .87, 'b', fontsize=14, weight='extra bold', transform=ax2.transAxes)
# Plot Simulation Trajectories
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
ax2 = plt.subplot(222)
plot_sim_results(ax2, ang_sim_dat, 'time1', 'tang_l*e1')
plot_sim_fits(ax2, ang_sim_dat, sim_lin_reg_params_glass, 'time1', msd=False)
# plt.plot(time, angle, 'o', label=labels[num], color=markercolors[num], markersize=2, markeredgewidth=0)
# plt.plot(time, fit_angle,'-', color='k', lw=.9, zorder=3)
#plt.xlabel('Time ($10^{-4}$s)')
plt.ylabel('Radians')
#plt.title('Intensity = I$_{o}$')
ax2.text(0.53, .95, r'Simulation $I_{0}$', transform=ax2.transAxes,
va='top', ha='center', fontsize=7.5, usetex=True)
plt.xlim(0,6.5)
y_ticks = np.array([0,1,2,3])*np.pi
y_label = [r"0", r"1$\pi$", r"2$\pi$", r"3$\pi$"]
plt.yticks(y_ticks,y_label)
plt.setp(ax2.get_xmajorticklabels(), visible=False)
ylim_2 = ax2.get_ylim()
ax2.set_ylim(ymax=1.0*ylim_2[1])
ax4 = plt.subplot(224)
plot_sim_results(ax4, ang_sim_dat, 'time1', 'tang_l*e2')
plot_sim_fits(ax4, ang_sim_dat, sim_lin_reg_params_plate, 'time1', msd=False)
plt.xlabel(r'Time ($\mathregular{10}^{\mathregular{-4}}$s)')
plt.ylabel('Radians')
#plt.title('Intensity = 4I$_{o}$')
ax4.text(0.53, .95, r'Simulation 4$I_{0}$', transform=ax4.transAxes,
va='top', ha='center', fontsize=7.5, usetex=True)
plt.xlim(0,6.5)
y_ticks = np.array([0,5,10])*np.pi
y_label = [r"0", r"5$\pi$", r"10$\pi$"]
plt.yticks(y_ticks,y_label)
ylim_4 = ax4.get_ylim()
ax4.set_ylim(ymax=1.0*ylim_4[1])
# ax3.text(0.03, .87, 'c', fontsize=14, weight='extra bold', transform=ax3.transAxes)
# ax4.text(0.03, .87, 'd', fontsize=14, weight='extra bold', transform=ax4.transAxes)
plt.tight_layout(pad=0.07)
#plt.subplots_adjust(hspace=0.07)
#plt.draw()
add_axes_label_inches(ax1, (0.059,0.059), '(a)', fontsize=10)
add_axes_label_inches(ax3, (0.059,0.059), '(b)', fontsize=10)
add_axes_label_inches(ax2, (0.059,0.059), '(c)', fontsize=10)
add_axes_label_inches(ax4, (0.059,0.059), '(d)', fontsize=10)
# ang_trajs.set_size_inches(3.25, 2.6)
# plt.tight_layout(pad=0.07)
plt.show()
matplotlib.rcParams['xtick.labelsize'] = 'medium'
matplotlib.rcParams['ytick.labelsize'] = 'medium'
In [26]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
ang_trajs.savefig(save_dir+"\Fig_3.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
ang_trajs.savefig(save_dir+"\Fig_3.pdf", dpi=800)
In [27]:
f = open('Fig_3.pkl', 'r')
notebook_info, frame_rate, sp_gl_plot, sp_pl_plot, sp_gl, sp_pl, TAMSD_gl_plot, TAMSD_pl_plot, lin_reg_params_glass, xf_gl, yf_gl, lin_reg_params_plate, xf_pl, yf_pl, TAMSD_gl_1_x, TAMSD_gl, TAMSD_gl_1_prams, TAMSD_pl_1_x, TAMSD_pl, TAMSD_pl_1_prams, frame_rate = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [28]:
'''Define functions to help calculate the anglular trajectory'''
def angular_traj(data_frame,x_cent,y_cent):
grp_tracks=data_frame.groupby('track id')
ang_trajs=[]
for name,grp in grp_tracks:
single_traj_angle=[0]
for pos in range(len(grp)-1):
angle_diff=calc_angle(grp['x pos'].iloc[pos+1],grp['y pos'].iloc[pos+1],x_cent,y_cent)-calc_angle(grp['x pos'].iloc[pos],grp['y pos'].iloc[pos],x_cent,y_cent)
if angle_diff>180:
angle_diff=angle_diff-360
if angle_diff<-180:
angle_diff=angle_diff+360
single_traj_angle.append(angle_diff)
ang_trajs.append(single_traj_angle)
return ang_trajs
In [29]:
fit_function=lambda x,*prams: prams[0]*x**2+prams[1]*x
In [30]:
TAMSD_gl_plot=TAMSD_gl[5:]
TAMSD_pl_plot=TAMSD_pl[5:]
ang_tamsd=plt.figure(figsize=[5,4])
# ang_trajs_and_tamsd=plt.figure(figsize=[10,8])
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
ax1=plt.subplot(221)
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
labels=['L=0','L=1','L=2','L=3','L=4','L=5',]
for k in range(len(TAMSD_gl_plot)):
lag_time = TAMSD_gl_1_x[k+5]
TAMSD_points = np.array(TAMSD_gl[k+5][0])*(np.pi/180)**2
TAMSD_fit = [fit_function(v,TAMSD_gl_1_prams[k+5].x[0],TAMSD_gl_1_prams[k+5].x[1])*(np.pi/180)**2 for v in TAMSD_gl_1_x[k+5]]
plt.scatter(lag_time, TAMSD_points,label=labels[k],color=markercolors[k], s=1.5, edgecolor='none')
plt.plot(lag_time, TAMSD_fit, '-', lw=.9, color='k')
#plt.legend(loc='upper left')
plt.legend(loc='upper left', labelspacing=0.0, frameon=False, prop={'size':7})
plt.yscale('log')
plt.xscale('log')
plt.ylim([10**-3,10**4])
plt.xlim([10**-2,10**1.3])
plt.xlabel('Lag Time (s)')
plt.ylabel(r'$\mathregular{\langle Radians^2\!\rangle}$')
#plt.setp(ax1.get_xticklabels(), visible=False)
plt.minorticks_off()
ax1.xaxis.labelpad=0
ax1.yaxis.labelpad=0
plt.title('Over Glass')
ax2=plt.subplot(222)
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
labels=['L=0','L=1','L=2','L=3','L=4','L=5',]
for k in range(len(TAMSD_pl_plot)):
lag_time = TAMSD_pl_1_x[k+5]
TAMSD_points = np.array(TAMSD_pl[k+5][0])*(np.pi/180)**2
TAMSD_fit = [fit_function(v,TAMSD_pl_1_prams[k+5].x[0],TAMSD_pl_1_prams[k+5].x[1])*(np.pi/180)**2 for v in TAMSD_pl_1_x[k+5]]
plt.scatter(lag_time, TAMSD_points,label=labels[k],color=markercolors[k], s=1.5, edgecolor='none')
plt.plot(lag_time, TAMSD_fit, '-', lw=.9, color='k')
plt.yscale('log')
plt.xscale('log')
plt.ylim([10**-3,10**5.5])
plt.xlim([10**-2,10**1.3])
plt.xlabel('Lag Time (s)')
plt.ylabel(r'$\mathregular{\langle Radians^2\!\rangle}$')
#ax2.text(-.1,1.1,r'$\bf{d}$', transform=ax2.transAxes)
ax2.xaxis.labelpad=0
ax2.yaxis.labelpad=0
plt.minorticks_off()
plt.title('Over Au Plate')
ax1.text(0.03, 0.87, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
ax2.text(0.03, 0.87, 'b', fontsize=14, weight='extra bold', transform=ax2.transAxes)
plt.tight_layout()
#plt.subplots_adjust(wspace= 0.4, hspace= 0.3)
#sns.despine()
plt.close()
plt.show()
Dynamic information about single Ag nanoparticles in the optical vortex trap for various substrates and angular drive l. Top panels show the trajectory of single Ag nanoparticle rotation with time over glass (a) and over the Au nanoplate (b). Bottom panels show the time averaged mean squared displacement of the single particle rotating over glass (c) and over the Au nanoplate (d). The dashed black curves in (a) and (b) show linear fits of the rotation versus time while the dashed black curves in (c) and (d) show fits to the TAMSD using eq (1).
Angular Speed and TAMSD Simulations¶
This is from Nishant's 2 particle simulations of the ring trap.
In [31]:
import scipy.io as sio
import re
ang_sim_dat = sio.loadmat('angular_displacement_msd.mat')
ang_msd_fit_dat = sio.loadmat('msd_fitting_params.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_sim_dat
new_ang_sim_dat = {}
for key, val in ang_sim_dat.iteritems():
if key in skip_list:
continue
new_ang_sim_dat[key] = val[:,0]
ang_sim_dat = new_ang_sim_dat
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_ang_msd_fit_dat = {}
for key, val in ang_msd_fit_dat.iteritems():
if key in skip_list:
continue
new_ang_msd_fit_dat[key] = val[0, :]
ang_msd_fit_dat = new_ang_msd_fit_dat
Fit the Simulation Data¶
In [32]:
import scipy
# Define fitting function for TAMSD simulations
fit_function=lambda x,*prams: prams[0]*x**2+prams[1]*x
def fit_function(x, a, b):
return a*x**2 + b*x
def curve_fit_positive(x_data,y_data):
def err_func((a,b)):
fit_y = a*x_data**2 + b*x_data
return ((fit_y-y_data)**2).sum()
a = scipy.optimize.curve_fit(fit_function, x_data, y_data, p0=(10**8, 50))
#a=scipy.optimize.minimize(err_func,(1,1),bounds=[(0,None),(0,None)],tol=0.0001 ,method='SLSQP')
return a[0]
# Fit TAMSD simulation over glass
sim_TAMSD_fit_glass = []
for L in range(6):
x = ang_sim_dat['time2'][:200]
y = ang_sim_dat['msd_l'+str(L)+'e1'][:200]
sim_TAMSD_fit_glass.append(curve_fit_positive(x, y))
# Fit TAMSD simulation over glass
sim_TAMSD_fit_plate = []
for L in range(6):
x = ang_sim_dat['time2'][:50]
y = ang_sim_dat['msd_l'+str(L)+'e2'][:50]
sim_TAMSD_fit_plate.append(curve_fit_positive(x, y))
In [33]:
import fnmatch
def plot_sim_results(axis, ang_sim_dat, time_key, regex):
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
keys = sorted(ang_sim_dat.keys())
match_keys = fnmatch.filter(keys, regex)
for num, key in enumerate(match_keys):
axis.plot(ang_sim_dat[time_key], ang_sim_dat[key], 'o', rasterized=True, c=markercolors[num], markersize=1.7, markeredgewidth=0)
axis.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
axis.xaxis.labelpad=0
axis.yaxis.labelpad=0
def plot_sim_fits(axis, ang_sim_dat, fit_params, time_key, msd=False):
x = ang_sim_dat[time_key]
if msd == False:
fit_function = lambda x,*prams: x*prams[0]
real_fit_params = fit_params
for L in range(6):
axis.plot(x, fit_function(x, real_fit_params[L]), '-k', rasterized=True)
elif msd == True:
fit_function=lambda x,*prams: prams[0]*x**2+prams[1]*x
real_fit_params = []
for v in fit_params:
real_fit_params.append(v)
for L in range(6):
axis.plot(x, fit_function(x, real_fit_params[L][0], real_fit_params[L][1]), '-k', rasterized=True, lw=0.7)
In [34]:
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
plt.figure(figsize=[5,4])
ax3 = plt.subplot(223)
plot_sim_results(ax3, ang_sim_dat, 'time2', 'msd_l*e1')
plot_sim_fits(ax3, ang_sim_dat, sim_TAMSD_fit_glass, 'time2', msd=True)
plt.xlabel('Lag Time (s)')
plt.ylabel(r'$\mathregular{\langle Radians^2\!\rangle}$')
plt.xlim(1.5*10**-7,2.5*10**-4)
plt.yscale('log')
plt.xscale('log')
plt.minorticks_off()
plt.title('Intensity = I$_{o}$')
ax4 = plt.subplot(224)
plot_sim_results(ax4, ang_sim_dat, 'time2', 'msd_l*e2')
plot_sim_fits(ax4, ang_sim_dat, sim_TAMSD_fit_plate, 'time2', msd=True)
plt.xlabel('Lag Time (s)')
plt.ylabel(r'$\mathregular{\langle Radians^2\!\rangle}$')
plt.xlim(1.5*10**-7,2.5*10**-4)
plt.yscale('log')
plt.xscale('log')
plt.minorticks_off()
plt.title('Intensity = 4I$_{o}$')
ax3.text(0.03, 0.87, 'c', fontsize=14, weight='extra bold', transform=ax3.transAxes)
ax4.text(0.03, 0.87, 'd', fontsize=14, weight='extra bold', transform=ax4.transAxes)
plt.tight_layout()
plt.close()
Angular Force vs L Experiment¶
In [35]:
f = open('Fig_4.pkl', 'r')
notebook_info, TAMSD_gl_1_prams, TAMSD_pl_1_prams, FL_vs_L_reg_params_pl, FL_vs_L_reg_params_gl, x = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [36]:
L_list=np.array([-5,-4,-3,-2,-1,0,1,2,3,4,5])
L_vs_ang_force=plt.figure(figsize=[3,2.5])
ax5 = plt.subplot(111)
ang_force_glass=np.sqrt(np.array([v.x[0] for v in TAMSD_gl_1_prams]))
ang_force_glass[:5]=ang_force_glass[:5]*-1
ang_force_plate=np.sqrt(np.array([v.x[0] for v in TAMSD_pl_1_prams]))
ang_force_plate[:5]=ang_force_plate[:5]*-1
ax5.scatter(L_list,ang_force_glass*(np.pi/180),s=20,marker='^', color="#4C72B0")
ax5.scatter(L_list,ang_force_plate*(np.pi/180), s=20, color="#C44E52", edgecolor='none')
ax5.plot(x,(x*FL_vs_L_reg_params_gl[0]+FL_vs_L_reg_params_gl[1])*(np.pi/180), color='k')
ax5.plot(x,(x*FL_vs_L_reg_params_pl[0]+FL_vs_L_reg_params_pl[1])*(np.pi/180), color='k')
#ax5.set_yticks([-2,-1,0,1,2])
plt.tick_params(axis='both', which='major')
plt.xlabel('Angular Number, L')
plt.ylabel(r'$\mathregular{\frac{{F_\mathit{l}}}{\gamma}}$ $\mathregular{(\frac{rad}{sec})}$')
plt.tight_layout()
plt.close()
Fit parameters of $\frac{F_l}{\gamma}$ from the TAMSD of single nanoparticles in the optical ring vortex over glass (blue) and over the gold nanoplate (red).
Angular Force vs L Simulations¶
From Nishant on this file:
File: msd_fitting_params In this file the variables fcex1 and fcex2 are the values of the fitting coefficients (fitting the msd to Eq 1 of the manuscript). The variable l denotes the different L values. So plot l vs fcex*. Also, x, y1 and x, y2 are linear fits to l, fcex1 and l, fcex2, respectively.
In [37]:
import scipy.io as sio
import re
ang_msd_fit_dat = sio.loadmat('msd_fitting_params.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_ang_msd_fit_dat = {}
for key, val in ang_msd_fit_dat.iteritems():
if key in skip_list:
continue
new_ang_msd_fit_dat[key] = val[0, :]
ang_msd_fit_dat = new_ang_msd_fit_dat
Used ratio of the slope of the over plate experiment result to determine the limits of the over plate simulation in order to get the graphs to line up.
In [38]:
ax6 = plt.subplot()
ax6.scatter(ang_msd_fit_dat['l'],ang_msd_fit_dat['fcex1']/10**3, marker='^', facecolor='none', color="#4C72B0")
ax6.scatter(ang_msd_fit_dat['l'],ang_msd_fit_dat['fcex2']/10**3, marker='o', facecolor='none', color="#C44E52", edgecolor="#C44E52")
ax6.plot(ang_msd_fit_dat['x'], ang_msd_fit_dat['y1']/10**3, color='#999999')
ax6.plot(ang_msd_fit_dat['x'], ang_msd_fit_dat['y2']/10**3, 'k')
ax6.set_xlabel('Angular Number, L')
ax6.set_ylabel(r'$\mathregular{\frac{{F_\mathit{l}}}{\gamma}}$ $\mathregular{(10^{3}\frac{rad}{sec})}$')
plt.close()
Both Same Plot¶
In [39]:
# Experimental Data
L_list=np.array([-5,-4,-3,-2,-1,0,1,2,3,4,5])
ang_force_glass=np.sqrt(np.array([v.x[0] for v in TAMSD_gl_1_prams]))
ang_force_glass[:5]=ang_force_glass[:5]*-1
ang_force_plate=np.sqrt(np.array([v.x[0] for v in TAMSD_pl_1_prams]))
ang_force_plate[:5]=ang_force_plate[:5]*-1
L_vs_ang_force=plt.figure(figsize=[3.5,2.5])
ax5 = plt.subplot(111)
ax5.scatter(L_list,ang_force_glass*(np.pi/180),s=20,marker='^', color="#4C72B0")
ax5.scatter(L_list,ang_force_plate*(np.pi/180), s=20, color="#C44E52", edgecolor='none')
ax5.plot(x,(x*FL_vs_L_reg_params_gl[0]+FL_vs_L_reg_params_gl[1])*(np.pi/180), color='k')
ax5.plot(x,(x*FL_vs_L_reg_params_pl[0]+FL_vs_L_reg_params_pl[1])*(np.pi/180), color='k')
#ax5.set_yticks([-2,-1,0,1,2])
ax5.set_ylim(-35,35)
plt.tick_params(axis='both', which='major')
plt.xlabel('Angular Number, L')
plt.ylabel(r'$\mathregular{\frac{{F_\mathit{l}}}{\gamma}}$ $\mathregular{(\frac{rad}{sec})}$')
plt.tight_layout()
# Simulation Data
ax6 = ax5.twinx()
ax6.scatter(ang_msd_fit_dat['l'],ang_msd_fit_dat['fcex1']/10**3, marker='^', facecolor='none', color="#4C72B0")
ax6.scatter(ang_msd_fit_dat['l'],ang_msd_fit_dat['fcex2']/10**3, marker='o', facecolor='none', color="#C44E52", edgecolor="#C44E52")
ax6.plot(ang_msd_fit_dat['x'], ang_msd_fit_dat['y1']/10**3, color='#999999')
ax6.plot(ang_msd_fit_dat['x'], ang_msd_fit_dat['y2']/10**3, 'k')
ax6.set_ylim(-68.5811,68.5811)
ax6.set_ylabel(r'$\mathregular{\frac{{F_\mathit{l}}}{\gamma}}$ $\mathregular{(10^{3}\frac{rad}{sec})}$')
plt.tight_layout()
plt.close()
Combine All Plots¶
In [40]:
from matplotlib.gridspec import GridSpec
matplotlib.rcParams['xtick.labelsize'] = 'small'
matplotlib.rcParams['ytick.labelsize'] = 'small'
TAMSD_gl_plot=TAMSD_gl[5:]
TAMSD_pl_plot=TAMSD_pl[5:]
fig_height = 2.8
fig_width = 7
ang_tamsd_and_force=plt.figure(figsize=[fig_width,fig_height])
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5',]
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
gs1 = GridSpec(2,2, hspace=0.1)
gs1.update(left=0.00, right=0.55)
# maintain_height_inches = 4.0
# maintain_width_inches = 5.0
# fig_coords_height = maintain_height_inches/fig_height
# fig_coords_width = maintain_width_inches/fig_width
# print fig_coords_height, fig_coords_width
# gs1 = GridSpec(2,2, wspace=0, hspace=0.35)
ax1=plt.subplot(gs1[0])
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
labels=['$l$=0','$l$=1','$l$=2','$l$=3','$l$=4','$l$=5']
for k in range(len(TAMSD_gl_plot)):
lag_time = TAMSD_gl_1_x[k+5]
TAMSD_points = np.array(TAMSD_gl[k+5][0])*(np.pi/180)**2
TAMSD_fit = [fit_function(v,TAMSD_gl_1_prams[k+5].x[0],TAMSD_gl_1_prams[k+5].x[1])*(np.pi/180)**2 for v in TAMSD_gl_1_x[k+5]]
plt.plot(lag_time, TAMSD_points,'o', rasterized=True,label=labels[k],color=markercolors[k], ms=1.7, markeredgewidth=0)
# plt.scatter(lag_time, TAMSD_points,label=labels[k],color=markercolors[k], s=1.5, edgecolor='none')
plt.plot(lag_time, TAMSD_fit, '-', rasterized=True, lw=.7, color='k')
leg_handles, leg_labels = ax1.get_legend_handles_labels()
leg_ax1 = plt.legend(leg_handles[::-1], leg_labels[::-1], loc='upper left',
labelspacing=-0.1, frameon=False, prop={'size':7},
bbox_to_anchor=(-0.1,.92), bbox_transform=ax1.transAxes,
handletextpad=-0.3)
leg_texts1 = leg_ax1.get_texts()
for text_obj in leg_texts1:
text_obj.set_usetex(True)
plt.yscale('log')
plt.xscale('log')
plt.ylim([10**-3,10**4])
plt.xlim([10**-2.2,10**1.3])
plt.yticks([10**-3, 10**-1, 10**1, 10**3], [r'$\mathregular{10}^{\mathregular{-3}}$', r'$\mathregular{10}^{\mathregular{-1}}$', r'$\mathregular{10}^{\mathregular{1}}$', r'$\mathregular{10}^{\mathregular{3}}$'])
#plt.xlabel('Lag Time (s)')
plt.ylabel(r'$\langle \mathregular{Radians}^\mathregular{2}\!\rangle}$')
plt.setp(ax1.get_xticklabels(), visible=False)
plt.minorticks_off()
ax1.xaxis.labelpad=0
ax1.yaxis.labelpad=0
#plt.title('Over Glass')
ax1.text(0.52, .94, 'Over Glass', transform=ax1.transAxes,
va='top', ha='center', fontsize=8, usetex=True)
# ax1.text(10**1.1, 10**-0.9, labels[0], transform=ax1.transData, va='center', ha='left', fontsize=7)
# ax1.text(10**1.1, 10**1.2, labels[1], transform=ax1.transData, va='center', ha='left', fontsize=7)
# ax1.text(10**1.1, 10**2.5, labels[2], transform=ax1.transData, va='center', ha='left', fontsize=7)
# ax1.text(10**1.1, 10**3, labels[3], transform=ax1.transData, va='center', ha='left', fontsize=7)
# ax1.text(10**1.1, 10**3.4, labels[4], transform=ax1.transData, va='center', ha='left', fontsize=7)
# ax1.text(10**1.1, 10**3.7, labels[5], transform=ax1.transData, va='center', ha='left', fontsize=7)
# xlim_1 = ax1.get_xlim()
# ax1.set_xlim(xmax=1.8*xlim_1[1])
ax2=plt.subplot(gs1[2])
markercolors=["#4C72B0", "#55A868", "#C44E52",
"#8172B2", "#CCB974", "#64B5CD"]
labels=['L=0','L=1','L=2','L=3','L=4','L=5',]
for k in range(len(TAMSD_pl_plot)):
lag_time = TAMSD_pl_1_x[k+5]
TAMSD_points = np.array(TAMSD_pl[k+5][0])*(np.pi/180)**2
TAMSD_fit = [fit_function(v,TAMSD_pl_1_prams[k+5].x[0],TAMSD_pl_1_prams[k+5].x[1])*(np.pi/180)**2 for v in TAMSD_pl_1_x[k+5]]
plt.plot(lag_time, TAMSD_points,'o', rasterized=True,label=labels[k],color=markercolors[k], ms=1.7, markeredgewidth=0)
plt.plot(lag_time, TAMSD_fit, '-', rasterized=True, lw=.7, color='k')
plt.yscale('log')
plt.xscale('log')
plt.ylim([10**-3.2,10**5.5])
plt.xlim([10**-2.2,10**1.3])
plt.yticks([10**-3, 10**-1, 10**1, 10**3, 10**5], [r'$\mathregular{10}^{\mathregular{-3}}$', r'$\mathregular{10}^{\mathregular{-1}}$', r'$\mathregular{10}^{\mathregular{1}}$', r'$\mathregular{10}^{\mathregular{3}}$', r'$\mathregular{10}^{\mathregular{5}}$'])
plt.xticks([10**-2, 10**-1, 10**0, 10**1], [r'$\mathregular{10}^{\mathregular{-2}}$', r'$\mathregular{10}^{\mathregular{-1}}$', r'$\mathregular{10}^{\mathregular{0}}$', r'$\mathregular{10}^{\mathregular{1}}$'])
plt.xlabel(r'Lag Time, $\tau$ (s)', usetex=True)
plt.ylabel(r'$\langle \mathregular{Radians}^\mathregular{2}\!\rangle}$')
#ax2.text(-.1,1.1,r'$\bf{d}$', transform=ax2.transAxes)
ax2.xaxis.labelpad=0
ax2.yaxis.labelpad=0
plt.minorticks_off()
#plt.title('Over Au Plate')
ax2.text(0.52, .94, 'Over Au Plate', transform=ax2.transAxes,
va='top', ha='center', fontsize=8, usetex=True)
ylim_2 = ax2.get_ylim()
ax2.set_ylim(ymax=3*ylim_2[1])
# ax1.text(0.03, 0.87, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
# ax2.text(0.03, 0.87, 'b', fontsize=14, weight='extra bold', transform=ax2.transAxes)
ax3 = plt.subplot(gs1[1])
plot_sim_results(ax3, ang_sim_dat, 'time2', 'msd_l*e1')
plot_sim_fits(ax3, ang_sim_dat, sim_TAMSD_fit_glass, 'time2', msd=True)
#plt.xlabel('Lag Time (s)')
plt.ylabel(r'$\langle \mathregular{Radians}^\mathregular{2}\!\rangle}$')
plt.xlim(1*10**-7,2.5*10**-4)
plt.yscale('log')
plt.xscale('log')
plt.minorticks_off()
plt.setp(ax3.get_xticklabels(), visible=False)
#plt.title('Intensity = I$_{o}$')
ax3.text(0.52, .94, 'Simulation $I_{0}$', transform=ax3.transAxes,
va='top', ha='center', fontsize=8, usetex=True)
y_ticks_ax3 = ax3.get_yticks()[1:-1]
print y_ticks_ax3
plt.yticks([10**(np.log10(v)) for v in y_ticks_ax3],
[r'$\mathregular{10}^{\mathregular{'+str(int(np.log10(c)))+r'}}$' for c in y_ticks_ax3])
ax4 = plt.subplot(gs1[3])
plot_sim_results(ax4, ang_sim_dat, 'time2', 'msd_l*e2')
plot_sim_fits(ax4, ang_sim_dat, sim_TAMSD_fit_plate, 'time2', msd=True)
plt.xlabel(r'Lag Time, $\tau$ (s)', usetex=True)
plt.ylabel(r'$\langle \mathregular{Radians}^\mathregular{2}\!\rangle}$')
plt.xlim(1*10**-7,2.5*10**-4)
plt.yscale('log')
plt.xscale('log')
plt.xticks([10**-7, 10**-6, 10**-5, 10**-4], [r'$\mathregular{10}^{\mathregular{-7}}$', r'$\mathregular{10}^{\mathregular{-6}}$', r'$\mathregular{10}^{\mathregular{-5}}$', r'$\mathregular{10}^{\mathregular{-4}}$'])
#plt.title('Intensity = 4I$_{o}$')
label4 = ax4.text(0.52, .94, r'Simulation 4$I_{0}$', transform=ax4.transAxes,
va='top', ha='center', fontsize=8, usetex=True)
ylim_4 = ax4.get_ylim()
ax4.set_ylim(ymax=3*ylim_4[1])
y_ticks_ax4 = ax4.get_yticks()[1:-2]
print y_ticks_ax4
plt.yticks([10**(np.log10(v)) for v in y_ticks_ax4],
[r'$\mathregular{10}^{\mathregular{'+str(int(np.log10(c)))+r'}}$' for c in y_ticks_ax4])
plt.minorticks_off()
# ax3.text(0.03, 0.87, 'c', fontsize=14, weight='extra bold', transform=ax3.transAxes)
# ax4.text(0.03, 0.87, 'd', fontsize=14, weight='extra bold', transform=ax4.transAxes)
matplotlib.rcParams['xtick.labelsize'] = 'medium'
matplotlib.rcParams['ytick.labelsize'] = 'medium'
# Experimental Data
gs2 = GridSpec(1,1, wspace=0, hspace=0.35)
gs2.update(left=0.65, right=0.95)
L_list=np.array([-5,-4,-3,-2,-1,0,1,2,3,4,5])
ang_force_glass=np.sqrt(np.array([v.x[0] for v in TAMSD_gl_1_prams]))
ang_force_glass[:5]=ang_force_glass[:5]*-1
ang_force_plate=np.sqrt(np.array([v.x[0] for v in TAMSD_pl_1_prams]))
ang_force_plate[:5]=ang_force_plate[:5]*-1
ax5 = plt.subplot(gs2[0])
ax5.scatter(L_list,ang_force_glass*(np.pi/180),s=20,marker='^', color="#4C72B0", zorder=1)
ax5.scatter(L_list,ang_force_plate*(np.pi/180), s=20, color="#C44E52", zorder=1)
ax5.plot(x,(x*FL_vs_L_reg_params_gl[0]+FL_vs_L_reg_params_gl[1])*(np.pi/180), color='k', zorder=20)
ax5.plot(x,(x*FL_vs_L_reg_params_pl[0]+FL_vs_L_reg_params_pl[1])*(np.pi/180), color='k', zorder=20)
#ax5.set_yticks([-2,-1,0,1,2])
ax5.set_ylim(-35,35)
ax5.yaxis.labelpad=0
plt.tick_params(axis='both', which='major')
plt.xlabel('Angular Number, $l$', usetex=True)
plt.ylabel(r'$\mathit{F(l)}/\gamma$ $\left(\mathsf{\frac{rad}{sec}}\right)$', usetex=True)
# ax5.text(0.03, 0.92, 'e', fontsize=14, weight='extra bold', transform=ax5.transAxes)
# Add arrows to show the axis that corresponds to experimental data
ax5.annotate('', xy=(-4,-4.3633), xycoords='data',
xytext=(-6,-11.5), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='<|-',
connectionstyle='arc3,rad=0.4'))
ax5.annotate('', xy=(-4.1,-22.6893), xycoords='data',
xytext=(-6,-18.5), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='<|-',
connectionstyle='arc3,rad=-0.4'))
# Simulation Data
ax6 = ax5.twinx()
ax6.scatter(ang_msd_fit_dat['l'],(ang_msd_fit_dat['fcex1']/10**3), marker='^', facecolor='none', color="#4C72B0", zorder=1)
ax6.scatter(ang_msd_fit_dat['l'],(ang_msd_fit_dat['fcex2']/10**3), marker='o', facecolor='none', color="#C44E52", edgecolor="#C44E52", zorder=1)
ax6.plot(ang_msd_fit_dat['x'], (ang_msd_fit_dat['y1']/10**3), color='#999999', zorder=20)
ax6.plot(ang_msd_fit_dat['x'], (ang_msd_fit_dat['y2']/10**3), 'k', zorder=20)
ax6.set_ylim(-68.5811,68.5811)
ax6.set_ylabel(r'$F(l)/\gamma$ $\left(\mathsf{10}^{\mathsf{3}} \frac{\mathsf{rad}}{\mathsf{sec}}\right)$', usetex=True)
ax6.yaxis.labelpad=7
# Add arrows to show the axis that corresponds to simulation data
ax6.annotate('', xy=(4,14), xycoords='data',
xytext=(6,23), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='<|-',
connectionstyle='arc3,rad=0.4'))
ax6.annotate('', xy=(5,51), xycoords='data',
xytext=(6,44), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='<|-',
connectionstyle='arc3,rad=-0.4'))
gs1.tight_layout(ang_tamsd_and_force, rect=[0,0,0.53, .98], pad=0)
gs2.tight_layout(ang_tamsd_and_force, rect=[0.55,0,1,.95], pad=0)
gs1.update(hspace=0.10, wspace=0.51)
# Add text labels after tight layout
add_axes_label_inches(ax1, (0.068,0.068), '(a)', fontsize=10)
add_axes_label_inches(ax2, (0.068,0.068), '(b)', fontsize=10)
add_axes_label_inches(ax3, (0.068,0.068), '(c)', fontsize=10)
add_axes_label_inches(ax4, (0.068,0.068), '(d)', fontsize=10)
add_axes_label_inches(ax5, (0.068,0.068), '(e)', fontsize=10)
print ax5.get_position()
plt.show()
print FL_vs_L_reg_params_gl
print FL_vs_L_reg_params_pl
In [41]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
ang_tamsd_and_force.savefig(save_dir+"\Fig_4.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
ang_tamsd_and_force.savefig(save_dir+"\Fig_4.pdf", dpi=800)
Exerimental Velocity around Ring¶
Load data from Mov_11121401 the single particle ramping over the nanoplate (see Ana_16011101). Only the L=4 part of the trajectory is used because it looks the best
In [42]:
import cPickle
f = open("C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Data\Fig_7.pkl")
trajs_pl = cPickle.load(f)
f.close()
In [43]:
def avg_st_dev_velocity_vs_theta(data_frame, t, mag_slider=1.6, theta_bin_num=18, ax=None):
um_conv=6.5/60/mag_slider/2
copy = data_frame.drop_duplicates(subset=['frame','track id']).copy()
dist = lambda x: (((x.shift(-t)['x pos']-x.shift(t)['x pos'])**2 + (x.shift(-t)['y pos']-x.shift(t)['y pos'])**2)**(0.5))/2.0
grouped = copy.groupby(['track id'])
# Need to stack if only one group
if len(grouped) == 1:
result = copy.groupby(['track id']).apply(dist)
result = result.stack()
else:
result = copy.groupby(['track id']).apply(dist)
copy['displacement'] = result.reset_index(level='track id').drop('track id', axis=1)
copy = copy.dropna(axis=0, subset=['displacement']).copy()
copy['velocity'] = copy['displacement']*um_conv*90.0
theta_bins = np.linspace(0, 360, num=theta_bin_num+1)
theta_bin_mid_points = (np.roll(theta_bins,-1)[:-1] + theta_bins[:-1])/2.0
copy['theta_bin'] = pd.cut(copy['theta'], theta_bins)
avgs = copy[['velocity','theta_bin']].groupby('theta_bin').mean().values[:,0]
std_dev = copy[['velocity','theta_bin']].groupby('theta_bin').std().values[:,0]
avgs = np.roll(avgs, -5)
std_dev = np.roll(std_dev, -5)
ax.plot(theta_bin_mid_points/180.0, avgs, '-', color="#4C72B0", mec="#4C72B0")
ax.fill_between(theta_bin_mid_points/180.0, avgs+std_dev, avgs-std_dev, facecolor="#4C72B0", alpha=0.5, lw=0.0)
#plt.plot(theta_bin_mid_points, avgs-std_dev)
plt.ylabel(u'Velocity (µm/sec)')
x_ticks = np.array([0,0.5,1,1.5,2])
x_label = [r"$0$", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]
plt.xticks(x_ticks, x_label)
#plt.xlabel('Theta')
In [44]:
%matplotlib inline
ax = plt.subplot()
avg_st_dev_velocity_vs_theta(trajs_pl[9], 1, theta_bin_num=20, ax=ax)
plt.close()
plt.show()
Simulation of Velocity around Ring¶
velocity_position_L4.mat is the velocity of a single particle in an L=4 trap. From Nishant:
This file contains the velocity of a single particle in an L=4 trap. theta_l4 (radians) and vmag_l4 (m/s) are the raw velocity and position data. mean_vmag is the binned average and theta_bin is the bin location. std_dev is the standard deviation for the binned velocity data.
theta-bin is actually the bin mid points for the mean velocity values (mean_vmag)
In [45]:
'''Load the simulation .mat file'''
import scipy.io as sio
vel_sim_dat = sio.loadmat('velocity_position_L4.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_vel_sim_dat = {}
for key, val in vel_sim_dat.iteritems():
if key in skip_list:
continue
#print val.shape, key
if key in ['theta_bin']:
new_vel_sim_dat[key] = val[0, :]
else:
new_vel_sim_dat[key] = val[:, 0]
vel_sim_dat = new_vel_sim_dat
In [46]:
fig, ax = plt.subplots()
theta = vel_sim_dat['theta_bin']
mean_vel = vel_sim_dat['mean_vmag']*10**6
std_dev = vel_sim_dat['std_dev']*10**6
ax.plot(theta/np.pi, mean_vel, '-o', color="#4C72B0", mec="#4C72B0")
ax.fill_between(theta/np.pi, mean_vel+std_dev, mean_vel-std_dev, facecolor="#4C72B0", alpha=0.5, lw=0.0)
x_ticks = np.array([0,0.5,1,1.5,2])
x_label = [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]
plt.xticks(x_ticks, x_label)
plt.ylabel(u'Velocity (µm/sec)')
plt.xlabel('Theta')
plt.title("Simulations Over Plate L=4")
plt.close()
plt.show()
Import Simulation Force Map Data¶
The force map and electric field intensity data are in forcemap_l4.mat are for L=4 from the simulations. From Nishant:
This file contains the optical force vectors and positions, and the electric field intensity in the ring trap. force_x (force_y) is the x (y) component of the force. pos_x (pos_y) is the position along x (y) axis. field_intensity is the electric field intensity and x, y are just x and y axes for the field intensity. The idea here is to plot the field intensity as the background and make a quiver plot with the force vectors.
In [47]:
'''Load the simulation .mat file'''
import scipy.io as sio
force_sim_dat = sio.loadmat('forcemap_l4.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_force_sim_dat = {}
for key, val in force_sim_dat.iteritems():
if key in skip_list:
continue
#print val.shape, key
if key in ['x', 'y']:
new_force_sim_dat[key] = val[0, :]
else:
new_force_sim_dat[key] = val
force_sim_dat = new_force_sim_dat
In [48]:
def add_two_polar_coords(ax, aspect='equal'):
"""A helper function that adds 0pi and pi/2 axis labels on a twinx and
twin y axis.
:param ax: A mpl.axes object you want to add polar labels to
:param aspect: The aspect ratio you want the final plot to have
"""
ax_y_polar = ax.twinx()
ax_y_polar.set_ylim(ax.get_ylim())
ax_y_polar.set_yticks([0])
ax_y_polar.set_yticklabels(['$0\pi$'])
ax_x_polar = ax.twiny()
ax_x_polar.set_xlim(ax.get_xlim())
ax_x_polar.set_xticks([0])
ax_x_polar.set_xticklabels(['$\pi/2$'])
ax.set(adjustable='box-forced', aspect=aspect)
ax_y_polar.set(adjustable='box-forced', aspect=aspect)
ax_x_polar.set(adjustable='box-forced', aspect=aspect)
In [49]:
def add_four_polar_coords(ax, size=None):
"""A helper function that adds 0pi, pi/2, pi, and 3pi/2 axis labels
on the inside of an axis following the unit circle.
:param ax: A mpl.axes object you want to add polar labels to
"""
if size == None:
size = 'large'
# Add the polar labels
ax.text(0.955, 0.5, '$\mathbf{\mathregular{0}}$', fontdict=dict(color='white', size=size),
weight='extra bold', transform=ax.transAxes, va='center', ha='right')
ax.text(0.5, 0.97, '$\mathbf{\pi/\mathregular{2}}$', fontdict=dict(color='white', size=size),
weight='extra bold', transform=ax.transAxes, va='top', ha='center')
ax.text(0.045, 0.5, '$\mathbf{\pi}$', fontdict=dict(color='white', size=size),
weight='extra bold', transform=ax.transAxes, va='center', ha='left')
ax.text(0.5, 0.03, '$\mathbf{\mathregular{3}\pi/\mathregular{2}}$', fontdict=dict(color='white', size=size),
weight='extra bold', transform=ax.transAxes, va='bottom', ha='center')
# Change ticks in graph
xticklines = ax.get_xticklines()
yticklines = ax.get_yticklines()
xcenter = len(xticklines)/2
ycenter = len(yticklines)/2
xticklines[xcenter].set(color='white', markeredgewidth=2, markersize=5, zorder=0)
xticklines[xcenter-1].set(color='white', markeredgewidth=2, markersize=5)
yticklines[ycenter].set(color='white', markeredgewidth=2, markersize=5)
yticklines[ycenter-1].set(color='white', markeredgewidth=2, markersize=5)
In [50]:
%matplotlib inline
force_fig = plt.figure(figsize=(4,4))
# Plot electric field intensity
field_x = force_sim_dat['x']/10**3
field_y = force_sim_dat['y']/10**3
field_int = force_sim_dat['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
# Plot the force vectors
force_pos_x = force_sim_dat['pos_x']/10**3
force_pos_y = force_sim_dat['pos_y']/10**3
force_mag_x = force_sim_dat['force_x']/10**3
force_mag_y = force_sim_dat['force_y']/10**3
plt.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], color="#4D0E10")
#plt.quiver(force_pos_x[:, :], force_pos_y[:, :], force_mag_x[:, :], force_mag_y[:, :], color="#4D0E10")
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
add_two_polar_coords(plt.gca())
plt.close()
Import the Force Map data from the Simulation L=0¶
This data is from Nishant's simulation for L=0. From Nishant:
File: forcemap_l0.mat x, y, field_intensity variables are the same as above. pos_x, pos_y and force_x, force_y are the coordinates and components of the force vectors for L=0.
In [51]:
'''Load the simulation .mat file'''
import scipy.io as sio
force_sim_dat_L0 = sio.loadmat('forcemap_l0.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_force_sim_dat_L0 = {}
for key, val in force_sim_dat_L0.iteritems():
if key in skip_list:
continue
#print val.shape, key
if key in ['x', 'y']:
new_force_sim_dat_L0[key] = val[0, :]
else:
new_force_sim_dat_L0[key] = val
force_sim_dat_L0 = new_force_sim_dat_L0
In [52]:
# Plot the Force Map
force_fig = plt.figure(figsize=(4,4))
# Plot electric field intensity
field_x = force_sim_dat_L0['x']*10**6
field_y = force_sim_dat_L0['y']*10**6
field_int = force_sim_dat_L0['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
# Plot the force vectors
force_pos_x = force_sim_dat_L0['pos_x']*10**6
force_pos_y = force_sim_dat_L0['pos_y']*10**6
force_mag_x = force_sim_dat_L0['force_x']*10**6
force_mag_y = force_sim_dat_L0['force_y']*10**6
plt.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], color="#4D0E10")
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
add_two_polar_coords(plt.gca())
plt.close()
Combining the Plots¶
In [53]:
from matplotlib.gridspec import GridSpec
# matplotlib.rcParams['axes.labelsize'] = 12
# matplotlib.rcParams['xtick.labelsize'] = 'medium'
# matplotlib.rcParams['ytick.labelsize'] = 'medium'
vel_exp_and_sim_force = plt.figure(figsize=(7,2.33333))
gs = GridSpec(2,3, width_ratios=[1,1.5,1.5])
ax1 = plt.subplot(gs[0,0])
avg_st_dev_velocity_vs_theta(trajs_pl[9], 1, theta_bin_num=20, ax=ax1)
#plt.title("Experiment L=4")
plt.setp(ax1.get_xticklabels(), visible=False)
ax1.set_ylim(ymax=1.2*ax1.get_ylim()[1])
ax1.set_yticks([50, 100, 150, 200])
ax1.tick_params(axis='y', labelsize='small')
ax1.text(0.56, .92, 'Experiment $l$=4', transform=ax1.transAxes,
va='top', ha='center', fontsize=8, usetex=True)
ax1.set_ylabel(u'Speed \n(µm/sec)')
ax2 = plt.subplot(gs[1,0])
theta = vel_sim_dat['theta_bin']
mean_vel = vel_sim_dat['mean_vmag']*10**6
std_dev = vel_sim_dat['std_dev']*10**6
ax2.plot(theta/np.pi, mean_vel, '-', color="#4C72B0", mec="#4C72B0")
ax2.fill_between(theta/np.pi, mean_vel+std_dev, mean_vel-std_dev, facecolor="#4C72B0", alpha=0.5, lw=0.0)
x_ticks = np.array([0,0.5,1,1.5,2])
x_label = [r"0", r"$\pi$/2", r"$\pi$", r"3$\pi$/2", r"2$\pi$"]
plt.xticks(x_ticks, x_label)
plt.ylabel(u'Speed \n(µm/sec)')
plt.xlabel(r'$\theta$')
#plt.title("Simulation L=4")
ax2.set_ylim(ymax=1.2*ax2.get_ylim()[1])
ax2.tick_params(axis='y', labelsize='small')
ax2.set_yticks([2000, 3000, 4000, 5000])
ax2.text(0.56, .92, 'Simulation $l$=4', transform=ax2.transAxes,
va='top', ha='center', fontsize=8, usetex=True)
# Plot the Force Map L=0
ax3 = plt.subplot(gs[:,1])
# Plot electric field intensity
field_x = force_sim_dat_L0['x']*10**6
field_y = force_sim_dat_L0['y']*10**6
field_int = force_sim_dat_L0['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis', rasterized=True)
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
# Plot the force vectors
force_pos_x = force_sim_dat_L0['pos_x']*10**6
force_pos_y = force_sim_dat_L0['pos_y']*10**6
force_mag_x = force_sim_dat_L0['force_x']*10**6
force_mag_y = force_sim_dat_L0['force_y']*10**6
plt.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], rasterized=True, color="#4D0E10")
plt.xlabel(u'x (µm)')
plt.ylabel(u'y (µm)')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
ax3.yaxis.labelpad=0
#add_two_polar_coords(plt.gca())
ax3.set(adjustable='box-forced', aspect='equal')
add_four_polar_coords(ax3, size='medium')
# Plot Force Map L=4
ax4 = plt.subplot(gs[:,2])
# Plot electric field intensity
field_x = force_sim_dat['x']/10**3
field_y = force_sim_dat['y']/10**3
field_int = force_sim_dat['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis', rasterized=True)
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
# Plot the force vectors
force_pos_x = force_sim_dat['pos_x']/10**3
force_pos_y = force_sim_dat['pos_y']/10**3
force_mag_x = force_sim_dat['force_x']/10**3
force_mag_y = force_sim_dat['force_y']/10**3
ax4.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], rasterized=True, color="#4D0E10")
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])
plt.xlabel(u'x (µm)')
plt.ylabel(u'y (µm)')
ax4.yaxis.labelpad=0
#add_two_polar_coords(ax4)
ax4.set(adjustable='box-forced', aspect='equal')
add_four_polar_coords(ax4, size='medium')
plt.tight_layout()
gs.update(hspace=0.1)
add_axes_label_inches(ax1, (0.040,0.042), '(a)', fontsize=10)
add_axes_label_inches(ax2, (0.040,0.042), '(b)', fontsize=10)
text_c = add_axes_label_inches(ax3, (0.042,0.042), '(c)', fontsize=10, color='white')
text_d = add_axes_label_inches(ax4, (0.042,0.042), '(d)', fontsize=10, color='white')
ax3.annotate('', xy=(107,125), xycoords='axes points',
xytext=(24,0), textcoords='offset points',
arrowprops=dict(edgecolor='#C7C7C7', facecolor='#C7C7C7', arrowstyle='<|-|>',
connectionstyle='arc3,rad=0'))
ax3.annotate('E', xy=(0,0), xytext=(119,119), textcoords='axes points',
horizontalalignment='center', verticalalignment='center', color='#C7C7C7', weight='extra bold')
ax4.annotate('', xy=(107,125), xycoords='axes points',
xytext=(24,0), textcoords='offset points',
arrowprops=dict(edgecolor='#C7C7C7', facecolor='#C7C7C7', arrowstyle='<|-|>',
connectionstyle='arc3,rad=0'))
ax4.annotate('E', xy=(0,0), xytext=(119,119), textcoords='axes points',
horizontalalignment='center', verticalalignment='center', color='#C7C7C7', weight='extra bold')
#vel_exp_and_sim_force.set_size_inches(7, 2.333333)
gs.tight_layout(vel_exp_and_sim_force, pad=0.1)
gs.update(hspace=0.1)
gs.update(wspace=0.4)
In [54]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
vel_exp_and_sim_force.savefig(save_dir+"\Fig_5.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
vel_exp_and_sim_force.savefig(save_dir+"\Fig_5.pdf", dpi=800)
In [55]:
fig5_data = pd.read_csv("C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Data\Fig_5.csv", index_col=0)
This data is from:
Ana_15101501 - Analysis on all Optical Binding Data NN Distributions, Polar Plot of NN
In [56]:
# In order to calculate the midpoint between the 1st nearest neighbors
# I need to get the position of the particles for a given nn_id and track
# id indexed by their track id and nn_id. Find the first nearest neighbor's
# position and add it to the data frame.
df = fig5_data.copy()
theta_nn_ids = df[['frame', 'nn_id', 'track id']]
frame_particle_id_reindex = df.set_index(['frame', 'track id', 'nn_id'])
theta_values_nn = frame_particle_id_reindex.ix[[tuple(v) for v in theta_nn_ids.values]].loc[:,['x pos','y pos']].values
df['nn_x_val'] = theta_values_nn[:,0]
df['nn_y_val'] = theta_values_nn[:,1]
# Once the nearest neighbors positions are found the midpoints are added
df['nn_x_mid'] = (df['x pos'] + df.nn_x_val)/2.0
df['nn_y_mid'] = (df['y pos'] + df.nn_y_val)/2.0
# The angle at the midpoint is calculated in the same polar coordinates
# used in the rest of the data frame.
one_point = df.ix[0]
x_cent, y_cent = calc_cent_from_polar(one_point['x pos'], one_point['y pos'], one_point.r, one_point.theta)
df['angle_midpoint'] = calc_angle(df.nn_x_mid, df.nn_y_mid, x_cent, y_cent)
# Must filter out i<j to prevent double counting
valid_positions = df[df['track id'] < df['nn_id']]
# Use all particles (not just 1st nearest neighbord)
#valid_positions = valid_positions[valid_positions['nn_num'] == 1]
mag_slider = 1.6
um_conv = 6.5/60/2.0/mag_slider
hist, xedges, yedges = np.histogram2d(valid_positions['theta'], valid_positions['nn_dist']*um_conv,
bins=(75,35), range=[[0,360],[0,2.3]], normed=True)
theta, r = np.meshgrid(xedges,yedges, indexing='ij')
Plotting the data
In [57]:
def shiftedColorMap(cmap, start=0, midpoint=0.5, stop=1.0, name='shiftedcmap'):
'''
Function to offset the "center" of a colormap. Useful for
data with a negative min and positive max and you want the
middle of the colormap's dynamic range to be at zero
Input
-----
cmap : The matplotlib colormap to be altered
start : Offset from lowest point in the colormap's range.
Defaults to 0.0 (no lower ofset). Should be between
0.0 and 1.0.
midpoint : The new center of the colormap. Defaults to
0.5 (no shift). Should be between 0.0 and 1.0. In
general, this should be 1 - vmax/(vmax + abs(vmin))
For example if your data range from -15.0 to +5.0 and
you want the center of the colormap at 0.0, `midpoint`
should be set to 1 - 5/(5 + 15)) or 0.75
stop : Offset from highets point in the colormap's range.
Defaults to 1.0 (no upper ofset). Should be between
0.0 and 1.0.
'''
cdict = {
'red': [],
'green': [],
'blue': [],
'alpha': []
}
# regular index to compute the colors
reg_index = np.linspace(start, stop, 257)
# shifted index to match the data
shift_index = np.hstack([
np.linspace(0.0, midpoint, 128, endpoint=False),
np.linspace(midpoint, 1.0, 129, endpoint=True)
])
for ri, si in zip(reg_index, shift_index):
r, g, b, a = cmap(ri)
cdict['red'].append((si, r, r))
cdict['green'].append((si, g, g))
cdict['blue'].append((si, b, b))
cdict['alpha'].append((si, a, a))
newcmap = matplotlib.colors.LinearSegmentedColormap(name, cdict)
plt.register_cmap(cmap=newcmap)
return newcmap
In [58]:
def make_cbar_match_polar_axis_height(polar_ax, cbar_ax):
"""Makes a color bar that corresponds to a polar projection plot
the same height as the polar plot.
This function works by converting data coordinates to of the top
and bottom of the polar plot into figure coordinates. Once figure
coordinates are discovered for the polar plot the colorbar's height is
adjusted to match the polar plot. This function only works if your x
axis (theta axis) goes from 0 to 2pi. The key to get this to work is
forcing matplotlib to draw the polar plot once so that the data
coordinates can be properly converted into figure coordinates.
:param polar_ax: The mpl.axes object that contains the polar plot
:param cbar_ax: The mpl.axes object that contains the color bar
"""
plt.draw()
fig = polar_ax.get_figure()
ylim_max = polar_ax.get_ylim()[1]
theta_offset = ax.get_theta_offset()
x_90 = (np.pi/2) - theta_offset
top_point_disp = polar_ax.transData.transform_point([x_90, ylim_max])
top_point_fig = fig.transFigure.inverted().transform(top_point_disp)
bot_point_disp = polar_ax.transData.transform_point([x_90+np.pi, ylim_max])
bot_point_fig = fig.transFigure.inverted().transform(bot_point_disp)
cbar_pos = cbar_ax.ax.get_position()
cbar_ax.ax.set_position([cbar_pos.x0, bot_point_fig[1],
cbar_pos.x1 - cbar_pos.x0,
top_point_fig[1] - bot_point_fig[1]])
In [59]:
new_cmap = shiftedColorMap(matplotlib.cm.afmhot_r, midpoint=0.30, name='shifted')
fig = plt.figure(figsize=(4,4))
ax = plt.subplot(111, projection='polar')
im = ax.pcolormesh(np.pi*theta/180.0, r, hist, cmap=new_cmap, rasterized=True)
cbar = plt.colorbar(im, ticks=[0,0.005,0.010,0.015], use_gridspec=False, pad=0.07)#, shrink=0.30, aspect=8, panchor=(0.0,0.7))
#cbar.ax.set_yticklabels([0, 5, 10,15])
cbar.set_label('Probability Density')
ax.set_theta_offset(np.pi*90/180.0)
x_label = [r"$\frac{\pi}{8}$", r"$\frac{3\pi}{8}$", r"$\frac{5\pi}{8}$", r"$\frac{7\pi}{8}$",
r"$\frac{9\pi}{8}$", r"$\frac{11\pi}{8}$", r"$\frac{13\pi}{8}$", r"$\frac{15\pi}{8}$"]
xlocs = [(-3 + v*2) * np.pi/8.0 for v in range(8)]
ax.set_thetagrids(np.array(xlocs)*180/np.pi, x_label, frac=1.13)#, fontsize=20)
ylocs = [0.6,1.2,1.8]
ylabels = [str(v) for v in ylocs]
ax.set_ylim([0,2.3])
ax.set_rgrids(ylocs, ylabels)#, rpad=0.0), angle=65)
ax.xaxis.set_tick_params(pad=100)
plt.grid(alpha=0.5)
make_cbar_match_polar_axis_height(ax, cbar)
test = ax.get_xmajorticklabels()[0]
print test.get_label()
font_prop = test.get_font_properties()
font_prop.get_size_in_points()
plt.close()
Combine all plots¶
In [62]:
from matplotlib.gridspec import GridSpec
NN_intro_fig = plt.figure(figsize=(3.34646,2.485948797))
# Plot NN Particles w.r.t. theta
ax1 = plt.subplot(projection='polar')
new_cmap = shiftedColorMap(matplotlib.cm.afmhot_r, midpoint=0.30, name='shifted')
im = ax1.pcolormesh(np.pi*theta/180.0, r, hist, cmap=new_cmap, rasterized=True)
cbar = plt.colorbar(im, ticks=[0,0.005,0.010,0.015], use_gridspec=True, pad=0.15)#, shrink=0.30, aspect=8, panchor=(0.0,0.7))
#cbar.ax.set_yticklabels([0, 5, 10,15])
cbar.set_label('Probability Density')
ax1.set_theta_offset(np.pi*90/180.0)
x_label = [r"$\pi$/8", r"3$\pi$/8", r"5$\pi$/8", r"7$\pi$/8",
r"9$\pi$/8", r"11$\pi$/8", r"13$\pi$/8", r"15$\pi$/8"]
xlocs = [(-3 + v*2) * np.pi/8.0 for v in range(8)]
ax1.set_thetagrids(np.array(xlocs)*180/np.pi, x_label, frac=1.22)#, fontsize=20)
ylocs = [0.6,1.2,1.8]
ylabels = [str(v) for v in ylocs]
ax1.set_ylim([0,2.3])
ax1.set_rgrids(ylocs, ylabels)#, rpad=0.0), angle=65)
ax1.xaxis.set_tick_params(pad=100)
plt.grid(alpha=1)
thetagrids = ax1.get_xgridlines()
for line in thetagrids:
line.set_linewidth(1)
line.set_dashes([8,4])
radialgrids = ax1.get_ygridlines()
for line in radialgrids:
line.set_linewidth(1)
#line.set_dashes([3,4])
line.set_alpha(0.7)
x_frac_offset = ax1.get_xlim()[0]
x_frac_total = abs(ax1.get_xlim()[0])+ax1.get_xlim()[1]
print x_frac_offset*0.25
# The \perp and \parallel symbols are not aligned with the center and are
# offset by 0.015 to the right of the center in fractional axes coordinates.
ax1.text(.995, 0.5, r"$\perp$", transform=ax1.transAxes, fontsize=18, va='center', ha='left')
ax1.text(0.485, 1.03, r"$\parallel$", transform=ax1.transAxes, fontsize=18, va='bottom', ha='center')
ax1.text(-0.025, 0.5, r"$\perp$", transform=ax1.transAxes, fontsize=18, va='center', ha='right')
ax1.text(0.485, -0.03, r"$\parallel$", transform=ax1.transAxes, fontsize=18, va='top', ha='center')
# Add E-Field Polarization Label
ax1.annotate('', xy=(-10,17), xycoords='axes points',
xytext=(32,0), textcoords='offset points',
arrowprops=dict(edgecolor='red', facecolor='red', arrowstyle='<|-|>',
connectionstyle='arc3,rad=0', lw=2))
ax1.annotate('E', xy=(0,0), xytext=(6,8), textcoords='axes points',
horizontalalignment='center', verticalalignment='center', color='red', weight='extra bold', size=15)
plt.tight_layout(pad=0.08)
make_cbar_match_polar_axis_height(ax1, cbar)
In [63]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
NN_intro_fig.savefig(save_dir+"\Fig_6.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
NN_intro_fig.savefig(save_dir+"\Fig_6.pdf", dpi=800)
Nearest neighbor distance of particles plotted with respect to the particle location in the optical ring vortex. The radial coordinate is the nearest neighbor separation between two particles. The angular coordinate is the position around the ring trap where one particle is located. The color represents the probability density of seeing events in the radial and angular coordinates. This data is from a multi-particle experiment at l=2 over glass.
In [177]:
f = open('Fig_6.pkl', 'r')
notebook_info, data_gl_all_l, data_gl_low_l, data_gl_mid_l, data_gl_high_l, data_pl_all_l, data_pl_low_l, data_pl_mid_l, data_pl_high_l = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [184]:
def plot_par_and_per_nn_hist_just_data(par_and_perp_list, glass=True):
'''par_and_perp_list features the parallel data first and the perpendicular
data second'''
if glass==True:
colors=['#162234','#4C72B0']
colors=list(reversed(colors))
if glass==False:
colors=['#551A12','#C44E52']
colors=list(reversed(colors))
axes=plt.gca()
# n_par,b_par,p_par=np.histogram(par_and_perp_list[0] ,
# bins=100, range=(0,3), normed=True)
# n_perp,b_perp,p_perp=plt.hist(par_and_perp_list[1] ,
# bins=100,range=(0,3),normed=True)
hist_par, edges = np.histogram(par_and_perp_list[0], bins=100, range=(0,3), normed=True)
hist_perp, edges = np.histogram(par_and_perp_list[1], bins=100, range=(0,3), normed=True)
midpoints = (edges[:-1] + edges[1:])/2.0
axes.plot(midpoints, hist_par, color=colors[0], linewidth=2, label=r"$\parallel$")
axes.plot(midpoints, hist_perp, color=colors[1], linewidth=1, label=r"$\perp$")
plt.xlabel(u'Separation (µm)')
#plt.ylabel('Probability Density')
axes.xaxis.labelpad=0
plt.xlim(0,3)
plt.xticks([0,0.6,1.2,1.8,2.4,3], [0,0.6,1.2,1.8,2.4,3])
#axes.get_xaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
axes.get_yaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
In [185]:
import matplotlib as mpl
# matplotlib.rcParams['xtick.labelsize'] = 'small'
# matplotlib.rcParams['ytick.labelsize'] = 'small'
low_mid_high_glass_plate_nn=plt.figure(figsize=[7,3.2941])
# ax_label = low_mid_high_glass_plate_nn.add_subplot(111)
# ax_label.set_ylabel("Probability Density", labelpad=12)
# plt.setp(ax_label.get_xticklabels(), visible=False)
# plt.setp(ax_label.get_yticklabels(), visible=False)
# for v in ax_label.spines.values():
# v.set_visible(False)
# # Over Glass Legend Lines
# parallel_g = mpl.lines.Line2D([], [], color='#4C72B0',
# linewidth=3,label=r"$\parallel$")
# perpendicular_g = mpl.lines.Line2D([], [], color='#162234',label=r"$\perp$")
# # Over Plate Legend Lines
# parallel_p = mpl.lines.Line2D([], [], color='#C44E52',
# linewidth=3,label=r"$\parallel$")
# perpendicular_p = mpl.lines.Line2D([], [], color='#551A12',label=r"$\perp$")
# Glass low Ls
ax1=low_mid_high_glass_plate_nn.add_subplot(231)
plot_par_and_per_nn_hist_just_data(data_gl_low_l)
ax1.set_title('Small $l$', usetex=True)
plt.ylabel('Probability\nDensity')#, fontsize=12)
ax1_handles, ax1_labels = ax1.get_legend_handles_labels()
ax1.legend(reversed(ax1_handles), reversed(ax1_labels), loc=(0.45,0.45), labelspacing=0.1, frameon=False)
plt.xlabel('')
# Glass mid Ls
ax2=low_mid_high_glass_plate_nn.add_subplot(232)
plot_par_and_per_nn_hist_just_data(data_gl_mid_l)
ax2.set_title('Medium $l$', usetex=True)
plt.xlabel('')
# Glass high Ls
ax3=low_mid_high_glass_plate_nn.add_subplot(233)
plot_par_and_per_nn_hist_just_data(data_gl_high_l)
ax3.set_title('Large $l$', usetex=True)
plt.xlabel('')
# Plate low Ls
ax4=low_mid_high_glass_plate_nn.add_subplot(234)
plot_par_and_per_nn_hist_just_data(data_pl_low_l, glass=False)
plt.ylabel('Probability\nDensity')
ax4_handles, ax4_labels = ax4.get_legend_handles_labels()
ax4.legend(reversed(ax4_handles), reversed(ax4_labels), loc=(0.45,0.45), labelspacing=0.1, frameon=False)
# Plate mid Ls
ax5=low_mid_high_glass_plate_nn.add_subplot(235)
plot_par_and_per_nn_hist_just_data(data_pl_mid_l, glass=False)
# Plate high Ls
ax6=low_mid_high_glass_plate_nn.add_subplot(236, sharex=ax4)
plot_par_and_per_nn_hist_just_data(data_pl_high_l, glass=False)
# Play witht the y limits so that the peaks don't touch the top of the plotting
# area.
ax1.set_ylim(0,4.6)
ax2.set_ylim(0,3.2)
ax3.set_ylim(0,3.6)
ax4.set_ylim(0,6.4)
ax5.set_ylim(0,3.8)
ax6.set_ylim(0,2.7)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
add_axes_label_inches(ax1, (0.068, 0.068), '(a)', fontsize=10)
add_axes_label_inches(ax2, (0.068, 0.068), '(b)', fontsize=10)
add_axes_label_inches(ax3, (0.068, 0.068), '(c)', fontsize=10)
add_axes_label_inches(ax4, (0.068, 0.068), '(d)', fontsize=10)
add_axes_label_inches(ax5, (0.068, 0.068), '(e)', fontsize=10)
add_axes_label_inches(ax6, (0.068, 0.068), '(f)', fontsize=10)
plt.tight_layout(pad=0.01)
plt.subplots_adjust(wspace= 0.2, hspace= 0.07)
plt.show()
In [186]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
low_mid_high_glass_plate_nn.savefig(save_dir+"\Fig_7.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
low_mid_high_glass_plate_nn.savefig(save_dir+"\Fig_7.pdf", dpi=800)
Nearest neighbor distributions grouped into all, small (L=0-2), medium (L=3), and large (L=4-5) applied optical force for experiments over glass (blue) and over the gold nanoplate (red). In each panel there are distributions for particles that are parallel (thick curve) or perpendicular (thin curve) to the polarization. The parallel distribution aggregates particle separations from the 3π⁄8 - 5π⁄8 and 11π⁄8 - 13π⁄8 regions of the optical ring trap while the perpendicular distribution aggregates particle separations from the π⁄8 - 15π⁄8 and 7π⁄8 - 9π⁄8 angular regions of the optical ring trap. Curves are plotted through the midpoints of the distribution bins.
fig8ac.mat is the 2 nanoparticle separation probability and potential from simulation for different drives. From Nishant:
fig8ac.mat contains the pair separation probability density function and the potentials of mean force (around the first optical binding separation) for different drives. low drive corresponds to l=0,1, medium l=2,3, and high l=4,5.
In [124]:
nn_sim_dat_full = mplhf.load_matlab_to_dict('fig8ac.mat')
fig8bd.mat is for 2 nanoparticle simulations that contain the PDF and PMF of interaction but looks at perpendicular (equator) or parallel (pole). From Nishant:
fig8bd.mat also contains pdfs and pmfs but for two different angular positions around the ring, 0 degree (equator) where the polarization is perpendicular to the nanoparticles and 90 degree (pole) where the polarization is parallel to the nanoparticles.
In [125]:
nn_eq_pol_sim_dat_full = mplhf.load_matlab_to_dict('fig8bd.mat')
In [126]:
matplotlib.rcParams['xtick.labelsize'] = 'small'
matplotlib.rcParams['ytick.labelsize'] = 'small'
fig = plt.figure(figsize=(3.34646,2.78868226))
colors = ["#2CA02C", "#9467BD", "#E377C2"]
# Panel 1: PDF NN Separation low, med, high
ax1 = plt.subplot(221)
plt.plot(nn_sim_dat_full['slow']/float(10**3), nn_sim_dat_full['pdlow']*10**2, label='Small $l$', color=colors[0])
#plt.plot(nn_sim_dat_full['sl']*10**6, nn_sim_dat_full['gfitlow'], color=colors[0], zorder=1)
plt.plot(nn_sim_dat_full['smed']/float(10**3), nn_sim_dat_full['pdmed']*10**2, label='Medium $l$', color=colors[1])
#plt.plot(nn_sim_dat_full['sm']*10**6, nn_sim_dat_full['gfitmid'], color=colors[1], zorder=1)
plt.plot(nn_sim_dat_full['shigh']/float(10**3), nn_sim_dat_full['pdhigh']*10**2, label='Large $l$', color=colors[2])
#plt.plot(nn_sim_dat_full['sh']*10**6, nn_sim_dat_full['gfithigh'], color=colors[2], zorder=1)
legend_ax1 = plt.legend(loc=(0.25,0.22), frameon=False, fontsize='small', labelspacing=0, handlelength=1)
legend_texts1 = legend_ax1.get_texts()
for text_obj in legend_texts1:
text_obj.set_usetex(True)
plt.xticks([0.0, 0.6, 1.2, 1.8])
plt.xlim(0.3,2.1)
plt.xlabel(u'Separation (µm)')
plt.ylabel('$\mathregular{10}^{\mathregular{-2}}$ PD')
#ax1.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
ax1.set_yticks(np.arange(0, 2.0, 0.4))
ax1.yaxis.labelpad = 0
ax1.xaxis.labelpad = 0
ax1.yaxis.major.formatter._useMathText = True
# Panel 2: PDF NN Separation equator, pole
ax2 = plt.subplot(222)
plt.plot(nn_eq_pol_sim_dat_full['sperp']/float(10**3), nn_eq_pol_sim_dat_full['pdperp']*10**2, label='$\perp$', color=colors[0])
#plt.plot(nn_sim_dat_full['sl']*10**6, nn_sim_dat_full['gfitlow'], color=colors[0], zorder=1)
plt.plot(nn_eq_pol_sim_dat_full['spara']/float(10**3), nn_eq_pol_sim_dat_full['pdpara']*10**2, label='$\parallel$', color=colors[1])
#plt.plot(nn_sim_dat_full['sm']*10**6, nn_sim_dat_full['gfitmid'], color=colors[1], zorder=1)
plt.legend(loc=(0.1,0.48), frameon=False, fontsize='small', labelspacing=0, handlelength=1)
plt.xticks([0.45, 0.55, 0.65, 0.75])
#plt.xlim(0.3,2.1)
plt.ylim(0.0,2.0)
plt.xlabel(u'Separation (µm)')
plt.ylabel('$\mathregular{10}^{\mathregular{-2}}$ PD')
#ax3.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
ax2.yaxis.labelpad = 0
ax2.xaxis.labelpad = 0
ax2.yaxis.major.formatter._useMathText = True
# Panel 3: PMF NN Separation low, med, high
ax3 = plt.subplot(223)
plt.plot(nn_sim_dat_full['slow']/float(10**3), nn_sim_dat_full['pmflow'], label='Low $l$', color=colors[0])
#plt.plot(nn_sim_dat_full['sl']*10**6, nn_sim_dat_full['gfitlow'], color=colors[0], zorder=1)
plt.plot(nn_sim_dat_full['smed']/float(10**3), nn_sim_dat_full['pmfmed'], label='Medium $l$', color=colors[1])
#plt.plot(nn_sim_dat_full['sm']*10**6, nn_sim_dat_full['gfitmid'], color=colors[1], zorder=1)
plt.plot(nn_sim_dat_full['shigh']/float(10**3), nn_sim_dat_full['pmfhigh'], label='High $l$', color=colors[2])
#plt.plot(nn_sim_dat_full['sh']*10**6, nn_sim_dat_full['gfithigh'], color=colors[2], zorder=1)
#plt.legend(loc='best', frameon=False)
plt.xlim(0.3,1.5)
plt.ylim(-0.25, 5.1)
ax3.set_xticks([0.3, 0.6, 0.9, 1.2, 1.5])
plt.xlabel(u'Separation (µm)')
plt.ylabel('PMF (units of $\mathregular{k}_\mathregular{B}\mathregular{T}$)')
ax3.yaxis.labelpad = 0
ax3.xaxis.labelpad = 0
# Panel 4: PMF NN Separation equator, pole
ax4 = plt.subplot(224)
plt.plot(nn_eq_pol_sim_dat_full['sperp']/float(10**3), nn_eq_pol_sim_dat_full['pmfperp'], label='$\perp$', color=colors[0])
#plt.plot(nn_sim_dat_full['sl']*10**6, nn_sim_dat_full['gfitlow'], color=colors[0], zorder=1)
plt.plot(nn_eq_pol_sim_dat_full['spara']/float(10**3), nn_eq_pol_sim_dat_full['pmfpara'], label='$\parallel$', color=colors[1])
#plt.plot(nn_sim_dat_full['sm']*10**6, nn_sim_dat_full['gfitmid'], color=colors[1], zorder=1)
#plt.legend(loc='best', frameon=False)
#plt.xticks([0.3, 0.6, 0.9])
#plt.xlim(0.36,.84)
plt.ylim(-0.2, 4.5)
plt.xticks([0.45, 0.55, 0.65, 0.75])
plt.xlabel(u'Separation (µm)')
plt.ylabel('PMF (units of $\mathregular{k}_\mathregular{B}\mathregular{T}$)')
ax4.yaxis.labelpad = 0
ax4.xaxis.labelpad = 0
mplhf.add_axes_label_inches(ax1, (0.055,0.055), '(a)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.055,0.055), '(c)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.055,0.055), '(b)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.055,0.055), '(d)', corner='upper right', fontsize=10)
plt.tight_layout(pad=0.00, h_pad=0.4, w_pad=0.4, rect=[0,0,1,0.995])#, h_pad=0.6)
matplotlib.rcParams['xtick.labelsize'] = 'medium'
matplotlib.rcParams['ytick.labelsize'] = 'medium'
print fig.get_size_inches()
In [127]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
fig.savefig(save_dir+"\Fig_8.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
fig.savefig(save_dir+"\Fig_8.pdf", dpi=800)
Figure 9¶
1D Langevin Simulations of Particle Pairs with Drive and Modulated Potentials¶
Import Nishant's Data¶
1DLDsimulation_modpairpot_drive.mat contains the data for the 1D Langevin dynamics simulation where particles are moved with different drives and have theta dependent angular drive that is either fixed or modulated. Inter particle potential is modulated in each.
1DLDsimulation_modpairpot_drive.mat contains the pdfs for different drives (l=0 to 5, and l=50) from a 1D Langevin simulation of two particles. scl0-scl50 (pdfscl0-pdfscl50) are the separations (pdfs) where the angular drive is assumed to be constant, while sml0-sml50 (pdfsml0-pdfsml50) are the separations (pdfs) assuming a modulated angular drive. For both cases, the inter particle potential is modulated.
In [136]:
'''Load simulation .mat file'''
import scipy.io as sio
langevin_1D_sim_dat_full = sio.loadmat('1DLDsimulation_modpairpot_drive.mat')
# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']
# Turn all the arrays into 1D arrays for the ang_sim_dat
new_langevin_1D_sim_dat_full = {}
for key, val in langevin_1D_sim_dat_full.iteritems():
if key in skip_list:
continue
new_langevin_1D_sim_dat_full[key] = val[0,:]
langevin_1D_sim_dat_full = new_langevin_1D_sim_dat_full
In [137]:
# Separate out the keys for the different types (constant, modulated, etc.)
langevin_1D_keys = langevin_1D_sim_dat_full.keys()
sep_const = sorted([key for key in langevin_1D_keys if 'scl' in key])[7:]
sep_mod = sorted([key for key in langevin_1D_keys if 'sml' in key])[7:]
pdf_const = sorted([key for key in langevin_1D_keys if 'pdfscl' in key])
pdf_mod = sorted([key for key in langevin_1D_keys if 'pdfsml' in key])
color_palette = ["#AEC7E8", "#CCB974", "#8172B2",
"#C44E52","#55A868","#4C72B0", 'k']
langevin_sim_fig = plt.figure(figsize=(3.34646,3.9661748))
ax1 = plt.subplot(211)
for num, (sep_key, pdf_key) in enumerate(zip(sep_const, pdf_const)):
if len(sep_key) == 5:
l_label = sep_key[-2:]
else:
l_label = sep_key[-1]
x = langevin_1D_sim_dat_full[sep_key]
y = langevin_1D_sim_dat_full[pdf_key]
plt.plot(x, y, color=color_palette[num], label=r'$l$='+l_label)
#ax1.set_yscale('log')
plt.xlabel('Separation (nm)')
plt.ylabel(r'$\mathregular{10}^\mathregular{-2}$ PD')
legend_ax1 = plt.legend(labelspacing=0.0, frameon=False, prop={'size':8.5})
legend_texts1 = legend_ax1.get_texts()
for text_obj in legend_texts1:
text_obj.set_usetex(True)
ax1.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
ax1.set_yticks(np.arange(0, 1.9, 0.3)*10**(-2))
#ax1.yaxis.major.formatter._useMathText = True
ax2 = plt.subplot(212)
for num, (sep_key, pdf_key) in enumerate(zip(sep_mod, pdf_mod)):
x = langevin_1D_sim_dat_full[sep_key]
y = langevin_1D_sim_dat_full[pdf_key]
plt.plot(x, y, color=color_palette[num])
#ax2.set_yscale('log')
plt.xlabel('Separation (nm)')
plt.ylabel(r'$\mathregular{10}^\mathregular{-2}$ PD')
ax2.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
ax2.set_yticks(np.arange(0, 1.9, 0.3)*10**(-2))
#ax2.yaxis.major.formatter._useMathText = True
add_axes_label_inches(ax1, (0.075, 0.074), '(a)', fontsize=10)
add_axes_label_inches(ax2, (0.075, 0.074), '(b)', fontsize=10)
plt.tight_layout(pad=0.1)
ax1.yaxis.get_offset_text().set_visible(False)
ax2.yaxis.get_offset_text().set_visible(False)
#ax1.text(0, 1.02, r'$\times\mathregular{10}^\mathregular{-2}$', transform=ax1.transAxes)
#ax2.text(0, 1.02, r'$\times\mathregular{10}^\mathregular{-2}$', transform=ax2.transAxes)
plt.show()
In [138]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
langevin_sim_fig.savefig(save_dir+"\Fig_9.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
langevin_sim_fig.savefig(save_dir+"\Fig_9.pdf", dpi=800)
In [63]:
f = open('Fig_6.pkl', 'r')
notebook_info, data_gl_all_l, data_gl_low_l, data_gl_mid_l, data_gl_high_l, data_pl_all_l, data_pl_low_l, data_pl_mid_l, data_pl_high_l = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [64]:
def plot_par_and_per_nn_hist_return_hist(par_and_perp_list, glass=True):
'''par_and_perp_list features the parallel data first and the perpendicular
data second and returns the midpoints, and two types of histograms.
:returns: midpoints (x), histogram parallel, histogram perpendicular
:rtype: arr, arr, arr
'''
if glass==True:
colors=['#162234','#4C72B0']
colors=list(reversed(colors))
if glass==False:
colors=['#551A12','#C44E52']
colors=list(reversed(colors))
axes=plt.gca()
# n_par,b_par,p_par=np.histogram(par_and_perp_list[0] ,
# bins=100, range=(0,3), normed=True)
# n_perp,b_perp,p_perp=plt.hist(par_and_perp_list[1] ,
# bins=100,range=(0,3),normed=True)
hist_par, edges = np.histogram(par_and_perp_list[0], bins=100, range=(0,3), normed=True)
hist_perp, edges = np.histogram(par_and_perp_list[1], bins=100, range=(0,3), normed=True)
midpoints = (edges[:-1] + edges[1:])/2.0
axes.plot(midpoints, hist_par, color=colors[0], linewidth=2, label=r"$\parallel$")
axes.plot(midpoints, hist_perp, color=colors[1], linewidth=1, label=r"$\perp$")
plt.xlabel(u'Separation (µm)')
#plt.ylabel('Probability Density')
axes.xaxis.labelpad=0
plt.xticks([0.6,1.2,1.8,2.4])
#axes.get_xaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
axes.get_yaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
return midpoints, hist_par, hist_perp
In [65]:
def add_inset_axes_coors(fig, axis, rect):
'''Returns a matplotlib axis object with the location set relative to an
input axis as a rectangle [left, up, width, height] in axis fraction
coordinates.
:params fig: The figure object you want to add the axis to
:params axis: The axis object (within fig) you want your new axis
coordinates relative to
:params rect: A rectangle defining your new axis in units of fraction of
input axis object of length 4. The parameters are [left, up, width, height].
'''
ax_x_start = rect[0]
ax_y_start = rect[1]
ax_x_end = rect[2] + ax_x_start
ax_y_end = rect[3] + ax_y_start
pix_xy_start = axis.transAxes.transform([ax_x_start, ax_y_start])
pix_xy_end = axis.transAxes.transform([ax_x_end, ax_y_end])
fig_width_height = fig.canvas.get_width_height()
rel_fig_start = pix_xy_start/fig_width_height
rel_fig_width_height = (pix_xy_end/fig_width_height) - rel_fig_start
inset_axis = fig.add_axes(np.concatenate((rel_fig_start,rel_fig_width_height)))
return inset_axis
In [66]:
# Glass low Ls
fig = plt.figure(figsize=(3.5,3))
ax1=plt.subplot()
midpoints_ag, hist_par_ag, hist_perp_ag = plot_par_and_per_nn_hist_return_hist(data_gl_low_l)
plt.legend([r"$\parallel$",r"$\perp$"],labelspacing=0.1, frameon=False)
plt.ylabel('Probability Density')
plt.xlabel(u'Separation (µm)')
# Plot Ag difference inset
ax1a = add_inset_axes_coors(fig, ax1, [0.5, 0.2, 0.45, 0.5])
ax1a.plot(midpoints_ag, hist_perp_ag-hist_par_ag, label="$\perp-\parallel$", color='#4c72b0')
ax1a.legend(frameon=False, labelspacing=0.1, fontsize=10, handletextpad=0, borderaxespad=0)
plt.ylim(-1, 3)
plt.xticks([0.6,1.2,1.8,2.4,3.0])
plt.yticks([-1,0,1.0,2.0])
plt.tight_layout()
plt.close()
Import the polystyrene over glass experiment¶
In [67]:
f = open('Fig_8.pkl', 'r')
notebook_info, par_and_perp_ps = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [68]:
import matplotlib as mpl
import scipy.ndimage
fig = plt.figure(figsize=(3.5,3))
['#162234','#4C72B0']
hist_par_ps, edges = np.histogram(par_and_perp_ps[0], bins=100, range=(0,3), normed=True)
hist_perp_ps, edges = np.histogram(par_and_perp_ps[1], bins=100, range=(0,3), normed=True)
midpoints_ps = (edges[:-1] + edges[1:])/2.0
ax2 = plt.subplot(111)
ax2.plot(midpoints_ps, hist_par_ps, '--', color='#4C72B0', label=r"$\parallel$", lw=2)
ax2.plot(midpoints_ps, hist_perp_ps, color='#162234', label=r"$\perp$", lw=1)
axes = plt.gca()
axes.get_xaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
axes.get_yaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
plt.legend(labelspacing=0.1, frameon=False)
plt.ylabel('Probability Density')
plt.xlabel(u'Separation (µm)')
# Plot PS difference inset
filtered = scipy.ndimage.filters.gaussian_filter1d(hist_perp_ps-hist_par_ps, 0.45)
ax2a = add_inset_axes_coors(fig, ax2, [0.5, 0.2, 0.45, 0.5])
ax2a.plot(midpoints_ps, filtered, label="$\perp-\parallel$", color='#4c72b0')
ax2a.legend(frameon=False, labelspacing=0.1, fontsize=10, handletextpad=0, borderaxespad=0)
plt.xticks([0.6,1.2,1.8,2.4,3.0])
plt.yticks([-0.3,0,0.3])
plt.tight_layout()
plt.close()
Combine both plots¶
In [70]:
import mpl_toolkits.axes_grid.inset_locator
ag_vs_glass_nn_sep = plt.figure(figsize=(3.34646,3.966175))
# PS L=1
hist_par_ps, edges = np.histogram(par_and_perp_ps[0], bins=100, range=(0,3), normed=True)
hist_perp_ps, edges = np.histogram(par_and_perp_ps[1], bins=100, range=(0,3), normed=True)
midpoints_ps = (edges[:-1] + edges[1:])/2.0
ax1 = plt.subplot(211)
ax1.plot(midpoints_ps, hist_perp_ps, color='#623c00', label=r"$\perp$", lw=1, zorder=1)
ax1.plot(midpoints_ps, hist_par_ps, '-', color='#F59700', label=r"$\parallel$", lw=2, zorder=0)
axes = plt.gca()
axes.get_xaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
axes.get_yaxis().set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
plt.legend(labelspacing=0.1, frameon=False, loc='lower left',
bbox_to_anchor=[0.17,0.45], bbox_transform=ax1.transAxes)
plt.ylabel('Probability Density')
plt.xlabel(u'Separation (µm)')
plt.xticks([0,0.6,1.2,1.8,2.4,3], [0,0.6,1.2,1.8,2.4,3])
plt.ylim(0,3.7)
# Ag Glass low Ls
ax2=plt.subplot(212)
midpoints_ag, hist_par_ag, hist_perp_ag = plot_par_and_per_nn_hist_return_hist(data_gl_low_l)
ax2_handles, ax2_labels = ax2.get_legend_handles_labels()
ax2.legend(reversed(ax2_handles), reversed(ax2_labels),
labelspacing=0.1, frameon=False,
loc='lower left', bbox_to_anchor=[0.17,0.45],
bbox_transform=ax2.transAxes)
plt.ylabel('Probability Density')
plt.xlabel(u'Separation (µm)')
plt.xticks([0,0.6,1.2,1.8,2.4,3], [0,0.6,1.2,1.8,2.4,3])
plt.tight_layout(pad=0.1, h_pad=1.1)
plt.draw()
# Plot PS difference inset
filtered = scipy.ndimage.filters.gaussian_filter1d(hist_perp_ps-hist_par_ps, 0.00)
ax1a = add_inset_axes_coors(ag_vs_glass_nn_sep, ax1, [0.52, 0.18, 0.44, 0.52])
ax1a.plot(midpoints_ps, filtered, label="$\perp-\parallel$", color='#8172B2', lw=1.5)
ax1a.legend(frameon=False, labelspacing=0.1, fontsize=10,
handletextpad=0, borderaxespad=0,
loc='lower left', bbox_to_anchor=[0.31,0.25],
bbox_transform=ax1a.transAxes)
ax1a.spines['top'].set_visible(False)
ax1a.spines['right'].set_visible(False)
ax1a.tick_params(axis='x', which='both', top='off')
ax1a.tick_params(axis='y', which='both', right='off')
ax1a.tick_params(labelsize='small')
plt.xticks([0,0.6,1.2,1.8,2.4,3], [0,0.6,1.2,1.8,2.4,3])
plt.yticks([-0.6,-0.3,0,0.3])
plt.yticks([-0.4,-0.2,0,0.2])
plt.ylim(-0.53,0.28)
# Plot Ag difference inset
ax2a = add_inset_axes_coors(ag_vs_glass_nn_sep, ax2, [0.52, 0.18, 0.44, 0.52])
ax2a.plot(midpoints_ag, hist_perp_ag-hist_par_ag, label="$\perp-\parallel$", color='#8172B2', lw=1.5)
ax2a.legend(frameon=False, labelspacing=0.1, fontsize=10,
handletextpad=0, borderaxespad=0,
loc='lower left', bbox_to_anchor=[0.27,0.4],
bbox_transform=ax2a.transAxes)
ax2a.spines['top'].set_visible(False)
ax2a.spines['right'].set_visible(False)
ax2a.tick_params(axis='x', which='both', top='off')
ax2a.tick_params(axis='y', which='both', right='off')
ax2a.tick_params(labelsize='small')
plt.ylim(-1, 2.5)
plt.xticks([0,0.6,1.2,1.8,2.4,3], [0,0.6,1.2,1.8,2.4,3])
plt.yticks([-1,0,1.0,2.0])
add_axes_label_inches(ax1, (0.075, 0.074), '(a)', fontsize=10)
add_axes_label_inches(ax2, (0.075, 0.074), '(b)', fontsize=10)
Out[70]:
In [71]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
ag_vs_glass_nn_sep.savefig(save_dir+"\Fig_10.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
ag_vs_glass_nn_sep.savefig(save_dir+"\Fig_10.pdf", dpi=800)
Nearest neighbor separations for 300nm polystyrene particles in the optical ring vortex over glass at l=1.
Supplemental Figures¶
Simulation Field of Ring over Au Plate¶
In [64]:
ring_field = mplhf.load_matlab_to_dict('fullnp_intensity_yz.mat')
In [65]:
ring_field.keys()
Out[65]:
In [66]:
efield_ring_plate = plt.figure(figsize=(3.34646, 2.2309733))
plt.pcolor(ring_field['y'], ring_field['z'], ring_field['Intensity_yz'].T, rasterized=True, cmap='viridis')
#high_int_args_z = np.argmax(ring_field['fyzsql1'].T[-70:,:137], axis=0)
#plt.plot(ring_field['y'][:137], ring_field['z'][-70+high_int_args_z], 'k', ls='--')
plt.ylabel('z (nm)')
plt.xlabel('y (nm)')
#plt.axis('equal')
print min(ring_field['y'])
plt.xlim(min(ring_field['y']), max(ring_field['y']))
plt.ylim(min(ring_field['z']), max(ring_field['z']))
plt.tight_layout(pad=0.01)
#plt.yticks(np.arange(200,-801, -200))
#plt.xticks(np.arange(-750, 801, 250))
In [67]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
efield_ring_plate.savefig(save_dir+"\Fig_S2.png", dpi=800)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\PhysRevE_OpticalRingVortex\main_text"
efield_ring_plate.savefig(save_dir+"\Fig_S2.pdf", dpi=800)
In [9]:
f = open('Fig_3.pkl', 'r')
notebook_info, frame_rate, sp_gl_plot, sp_pl_plot, sp_gl, sp_pl, TAMSD_gl_plot, TAMSD_pl_plot, lin_reg_params_glass, xf_gl, yf_gl, lin_reg_params_plate, xf_pl, yf_pl, TAMSD_gl_1_x, TAMSD_gl, TAMSD_gl_1_prams, TAMSD_pl_1_x, TAMSD_pl, TAMSD_pl_1_prams, frame_rate = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
Import simulation data for pmf of single particles around ring
In [10]:
pmf_sim = mplhf.load_matlab_to_dict('pmf_1np_ringtrap_lcmp.mat')
In [14]:
color_cycle = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd']
pmf_plate = plt.figure(figsize=(5,3))
ax1 = plt.subplot(121)
for idx, traj in enumerate(sp_pl[6:]):
L = idx+1
full_rotations = domf.find_full_trajs_around_ring(traj)
full_rotations = full_rotations[full_rotations['full_track_region']==True]
radians = ((full_rotations.theta - 90.0)%360) * np.pi/180.0
x_values, estimate = periodic_1D_guassian_kde(radians, bandwidth=0.1)
pmf = -np.log(estimate)
pmf = pmf - np.mean(pmf)
plt.plot(x_values, pmf, color_cycle[idx], label='$l$='+str(L), lw=1.5)
plt.title('Experiments Over Plate')
plt.ylabel('pmf ($\mathregular{k_B}$T)')
plt.xlabel('$\mathregular{\Theta}$ (radians)')
plt.xlim(0, 2*np.pi)
plt.legend(loc='best', labelspacing=-0.25, frameon=False, handlelength=1, fontsize='small', handletextpad=0.5, borderaxespad=0)
x_labels = [r"0", r"$\pi$/2", r"$\pi$", r"3$\pi$/2", r"2$\pi$"]
plt.xticks(np.linspace(0, 2*np.pi, 5), x_labels)
ax2 = plt.subplot(122)
for l in range(2,6):
pmf = pmf_sim['pmf'+str(l)+'s']
pmf = pmf - np.mean(pmf)
plt.plot(pmf_sim['t_values'], pmf, color_cycle[l-1], label='$l$='+str(l), lw=1.5)
plt.title('ED-LD Simulations 4$\mathregular{I_o}$')
plt.ylabel('pmf ($\mathregular{k_B}$T)')
plt.xlabel('$\mathregular{\Theta}$ (radians)')
plt.xlim(0, 2*np.pi)
#plt.legend(loc='best', labelspacing=0, frameon=False, handlelength=1, handletextpad=0.5)
x_labels = [r"0", r"$\pi$/2", r"$\pi$", r"3$\pi$/2", r"2$\pi$"]
plt.xticks(np.linspace(0, 2*np.pi, 5), x_labels)
mplhf.add_axes_label_inches(ax1, (0.080,0.05), '(a)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.080,0.05), '(b)', corner='upper right', fontsize=10)
plt.tight_layout(pad=0.01)
plt.show()
In [15]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
pmf_plate.savefig(save_dir+"\Fig_reply_1.png", dpi=600)
Extra Figures¶
Example NN Trajectory and Cumalitive Dwell Times for L=1 and L=4 Over Glass¶
Example NN Trajectory
In [158]:
f = open('Fig_7a.pkl', 'r')
notebook_info, time, nn_dist = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [159]:
f = open('Fig_7bc.pkl', 'r')
notebook_info, cum_dwell_L_1, cum_dwell_L_4 = cPickle.load(f)
f.close()
print "This data is from:\n"+notebook_info[0]
In [160]:
'''Plot the graphs and subplot here'''
first_state=[0.4,0.8]
second_state=[0.9,1.4]
third_state=[1.6,2.0]
colors = ["#55A868",
"#C44E52",
"#8172b2"]
fig = plt.figure(figsize=(7,5))
ax1 = plt.subplot2grid((3,2), (0,0), colspan=2, rowspan=1)
ax1.plot(time, nn_dist, color="#354f7b")
plt.ylim(0,2.6)
plt.yticks(np.arange(0, 3.0, 0.6))
#plt.ylim(0,3)
xmax = plt.xlim()[1]
plt.fill_between(np.arange(xmax+1),first_state[0],first_state[1], color=colors[0], alpha=0.3)
plt.fill_between(np.arange(xmax+1),second_state[0],second_state[1], color=colors[1],alpha=0.3)
plt.fill_between(np.arange(xmax+1),third_state[0],third_state[1], color=colors[2], alpha=0.5)
plt.xlabel('Time (s)')
plt.ylabel(u'Seperation (µm)')
#ax_label = plt.subplot2grid((2,3), (0,2), rowspan=2)
ax2 = plt.subplot2grid((3,2), (1,0), rowspan=2)
for num,i in enumerate(range(0,6,2)):
ax2.scatter(cum_dwell_L_1[i], cum_dwell_L_1[i+1], c=colors[num], edgecolor='none')
plt.yscale('log')
plt.ylim(10**-3, 2)
plt.xlim(0,0.45)
plt.xlim(0,0.45)
plt.xticks(np.arange(0, 0.5, 0.1))
plt.xlabel('Dwell Time (s)')
plt.ylabel('Cumulative Dwell Time')
#plt.setp(ax2.get_xticklabels(), visible=False)
plt.yticks([10**-2, 10**-1, 1])
ax3 = plt.subplot2grid((3,2), (1,1), rowspan=2, sharex=ax2)
for num,i in enumerate(range(0,6,2)):
#ax3.plot(cum_dwell_L_4[i], cum_dwell_L_4[i+1], 'o')
ax3.scatter(cum_dwell_L_4[i], cum_dwell_L_4[i+1], c=colors[num], edgecolor='none')
plt.yscale('log')
plt.ylim(10**-3, 2)
plt.xlim(0,0.45)
plt.xticks(np.arange(0, 0.5, 0.1))
plt.xlabel('Dwell Time (s)')
plt.yticks([10**-2, 10**-1, 1])
ax1.text(-0.13, 1.2, 'a', fontsize=14, weight='extra bold', transform=ax1.transAxes)
ax2.text(-0.15, 1.1, 'b', fontsize=14, weight='extra bold', transform=ax2.transAxes)
ax3.text(-0.15, 1.1, 'c', fontsize=14, weight='extra bold', transform=ax3.transAxes)
plt.tight_layout()
plt.subplots_adjust(hspace=1.2, wspace=0.2)
# matplotlib.rcParams
# print ax2.yaxis.get_label().get_position()
# sharey_x = ax2.yaxis.get_label().get_position()[0]
# ax2.yaxis.set_label_coords('left',0, transform=fig.transFigure)
In [161]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
fig.savefig(save_dir+"\Fig_7.png", dpi=200)
Analysis of nearest neighbor trajectories. (a) A typical trajectory of one particle moving relative to its nearest neighbor. The shaded regions show the first (green), second (red), and third (blue) optical binding sites. (b) The cumulative dwell time distribution of a nearest neighbor trajectory for small drive. (c) The cumulative dwell time distribution of a nearest neighbor trajectory for large drive. In (b) and (c) green, red, and blue, markers correspond to the first, second, and third optical binding regions respectfully.